多视角子空间学习系列之 MCCA (Multi-view CCA) 多视角CCA Horst算法

优化目标

前面已经讲了典型相关分析CCA，并且提到，CCA是一种双视角的方法，只能处理视角数为2的情况。为了将CCA应用于更多视角，一些研究人员提出了MCCA(Multi-view CCA)，即多视角CCA，将CCA直观地扩展到多视角版本。给定mmm个视角X1,X2,⋯,Xm{X_1,X_2,\cdots,X_m}X1,X2,⋯,Xm，其中第iii个视角Xi∈RDi×nX_i\in \mathbb{R}^{D_i\times n}Xi∈RDi×n，DiD_iDi为其维度，nnn为样本数量，各视角样本数量相同。MCCA希望找到W1,W2,⋯,Wm{W_1,W_2,\cdots,W_m}W1,W2,⋯,Wm将各视角投影到公共子空间，其中Wi∈RDi×dW_i\in \mathbb{R}^{D_i\times d}Wi∈RDi×d，ddd为子空间维度，使得子空间中任意两个视角之间的相关系数之和：
max⁡W1,⋯,Wm∑i,jmWiTXiXjTWjs.t.WiTXiXiTWi=I\max_{W_1,\cdots,W_m} \sum_{i,j}^m W_i^TX_iX_j^TW_j \\ s.t.\ W_i^TX_iX_i^T W_i=I W1,⋯,Wmmaxi,j∑mWiTXiXjTWjs.t. WiTXiXiTWi=I

本文仍然既是总结又是读论文笔记。

d=1d=1d=1的情况

d>1d>1d>1 的情况有点复杂，先看d=1d=1d=1的情况如何解。首先d=1d=1d=1时，为了方便，我们用w1,⋯,wm{w_1,\cdots,w_m}w1,⋯,wm代替W1,⋯,Wm{W_1,\cdots,W_m}W1,⋯,Wm，上式写为：
max⁡w1,⋯,wm∑i,jmwiTXiXjTwjs.t.wiTXiXiTwi=1(1)\max_{w_1,\cdots,w_m} \sum_{i,j}^m w_i^TX_iX_j^Tw_j\tag{1} \\ s.t.\ w_i^TX_iX_i^T w_i=1 w1,⋯,wmmaxi,j∑mwiTXiXjTwjs.t. wiTXiXiTwi=1(1)

令C(i,i)=XiXiTC^{(i,i)}=X_iX_i^TC(i,i)=XiXiT，对C(i,i)C^{(i,i)}C(i,i)做Cholesky分解得C(i,i)=DiTDiC^{(i,i)}=D_i^TD_iC(i,i)=DiTDi，其实Di=XiTD_i=X_i^TDi=XiT，然后令XiTwi=xiX_i^Tw_i=x_iXiTwi=xi，则wi=Di−1xiw_i=D_i^{-1}x_iwi=Di−1xi. 令：
Aij=Di−TC(i,j)Dj−1A_{ij}=D_i^{-T}C^{(i,j)}D_j^{-1} Aij=Di−TC(i,j)Dj−1

我就寻思这AijA_{ij}Aij其实是个单位矩阵。原文这么写有他的道理。往后看吧。然后公式(1)(1)(1)就可以写为：
max⁡x1,⋯,xm∑i,jmxiTAijxjs.t.xiTxi=1\max_{x_1,\cdots,x_m}\sum_{i,j}^m x_i^TA_{ij} x_j \\ s.t.\ x_i^Tx_i=1 x1,⋯,xmmaxi,j∑mxiTAijxjs.t. xiTxi=1

进一步写为：
max⁡xxTAxs.t.xiTxi=1\max_{x} x^TAx \\ s.t.\ x_i^Tx_i=1 xmaxxTAxs.t. xiTxi=1

这里的xxx是由x1,⋯,xmx_1,\cdots,x_mx1,⋯,xm这mmm个列向量拼成的列向量。但是算的时候又按照矩阵分块运算来来算的。我们来套一下Lagrangian乘子法：
L(x,λ)=∑i,jmxiTAijxj+∑imλi(1−xiTxi)L(x,\lambda)=\sum_{i,j}^m x_i^TA_{ij}x_j+\sum_i^m\lambda_i(1-x_i^Tx_i) L(x,λ)=i,j∑mxiTAijxj+i∑mλi(1−xiTxi)

对xix_ixi求导：
∂∂xiL(x,λ)=∑jmAijxj−2λixi=0∑jmAijxj=2λixi\frac{\partial}{\partial x_i}L(x,\lambda)=\sum_{j}^m A_{ij}x_j-2\lambda_ix_i=0 \\ \sum_{j}^m A_{ij}x_j=2\lambda_ix_i ∂xi∂L(x,λ)=j∑mAijxj−2λixi=0j∑mAijxj=2λixi

因此构成一个大的矩阵块的形式：
[A11⋯A1m⋮⋱⋮Am1⋯Amm][x1⋮xm]=[λ1x1⋮λmxm]\left[ \begin{array}{ccc} A_{11} & \cdots & A_{1m} \\ \vdots & \ddots & \vdots \\ A_{m1} & \cdots & A_{mm} \\ \end{array} \right] \left[ \begin{array}{c} x_{1} \\ \vdots \\ x_{m} \\ \end{array} \right]= \left[ \begin{array}{c} \lambda_1x_1 \\ \vdots \\ \lambda_mx_m \\ \end{array} \right] ⎣⎢⎡A11⋮Am1⋯⋱⋯A1m⋮Amm⎦⎥⎤⎣⎢⎡x1⋮xm⎦⎥⎤=⎣⎢⎡λ1x1⋮λmxm⎦⎥⎤

这里把所有2λi2\lambda_i2λi都写成λi\lambda_iλi的形式了，都是常数参数而已。注意本文对AAA矩阵的定义与原文不同。

然后该问题就变成了一个多元的广义特征值与特征向量问题。文章里面对于这种AAA为对称、正定的generic matrix（特征值都不相同的矩阵），使用了一种普遍的算法Horst算法：

步骤上是容易理解的，这是一种迭代算法，xxx需要有一个人为设置的初值，maxitermaxitermaxiter是人为设置的最大迭代次数。Horst算法的正确性和思路都有由来，但是我现在不想再细究了，论文里说这个算法能保证收敛到一种局部最优解。

形式转换

论文在这里做了一些形式转换，加了一些正则化考虑。公式(1)(1)(1)中，如果特征维度DiD_iDi小于样本数nnn，XiXiTX_iX_i^TXiXiT有可能是奇异的（不可逆的），而且有可能造成过拟合。因此文章将公式(1)(1)(1)改为如下形式：
max⁡w1,⋯,wm∑i,jmwiTXiXjTwjs.t.wiT(1−κDi−1XiXiT+κI)wi=1\max_{w_1,\cdots,w_m} \sum_{i,j}^m w_i^TX_iX_j^Tw_j \\ s.t.\ w_i^T (\frac{1-\kappa}{D_i-1}X_iX_i^T+\kappa I) w_i=1 w1,⋯,wmmaxi,j∑mwiTXiXjTwjs.t. wiT(Di−11−κXiXiT+κI)wi=1

其中κ∈[0,1]\kappa \in [0,1]κ∈[0,1]，然后令wi=Xiyi,Ki=XiTXiw_i=X_iy_i,K_i=X_i^TX_iwi=Xiyi,Ki=XiTXi，则yi=Xi−1wiy_i=X_i^{-1}w_iyi=Xi−1wi，然后上式等价于：
max⁡y∑i,jyiTKiKjTyjs.t.yiT(1−κDi−1KiKiT+κKi)yi=1\max_y \sum_{i,j} y_i^TK_iK_j^Ty_j \\ s.t.\ y_i^T (\frac{1-\kappa}{D_i-1}K_iK_i^T+\kappa K_i) y_i=1 ymaxi,j∑yiTKiKjTyjs.t. yiT(Di−11−κKiKiT+κKi)yi=1

推导上注意，矩阵转置的逆等于其逆的转置。文章认为这么做的好处是便于引入使用核函数。

进一步，文章认为KiK_iKi也是Typically病态的，甚至数据中心化后经常是奇异矩阵，并且处于限制方差和系数的原因，令：
Ki~=(1−κDi−1Ki+κ2Di−11−κI)\widetilde{K_i}=(\sqrt{\frac{1-\kappa}{D_i-1}}K_i+\frac{\kappa}{2}\sqrt{\frac{D_i-1}{1-\kappa}}I) Ki=(Di−11−κKi+2κ1−κDi−1I)

因此：
1−κDi−1KiKiT+κKi≈Ki~Ki~T\frac{1-\kappa}{D_i-1}K_iK_i^T+\kappa K_i \approx \widetilde{K_i}\widetilde{K_i}^T Di−11−κKiKiT+κKi≈KiKiT

因此最终把优化目标的形式变为：
max⁡y∑i,jyiTKiKjTyjs.t.yiTKi~Ki~Tyi=1\max_y \sum_{i,j} y_i^TK_iK_j^Ty_j \\ s.t.\ y_i^T \widetilde{K_i} \widetilde{K_i}^T y_i=1 ymaxi,j∑yiTKiKjTyjs.t. yiTKiKiTyi=1

d>1d>1d>1的情况

接上文。d>1d>1d>1时，优化目标为：
max⁡y∑i,jyiTKiKjTyjs.t.yiTKi~Ki~Tyi=1YiTKi~Ki~Tyi=0\max_y \sum_{i,j} y_i^TK_iK_j^Ty_j \\ s.t.\ y_i^T \widetilde{K_i} \widetilde{K_i}^T y_i=1 \\ Y_i^T \widetilde{K_i} \widetilde{K_i}^T y_i=0 ymaxi,j∑yiTKiKjTyjs.t. yiTKiKiTyi=1YiTKiKiTyi=0

这啥意思呢，yiy_iyi表示的是一个维度，因为d>1d>1d>1，子空间中各个维度应该是互不相关的，这里YiY_iYi表示的除了yiy_iyi之外的其他维度。

为了应用Horst算法，论文先做如下转换：
Zi=Ki~Yi,zi=Ki~yiZ_i=\widetilde{K_i} Y_i,z_i=\widetilde{K_i} y_i Zi=KiYi,zi=Kiyi

然后定义算子：
Pi=I−Ki~YiYiTKi~=I−ZiZiTP_i=I-\widetilde{K_i} Y_iY_i^T\widetilde{K_i} =I-Z_iZ_i^T Pi=I−KiYiYiTKi=I−ZiZiT

然后优化目标可以写为：
max⁡z∑i,jziTKi~−TKiKjTKj~−1zjs.t.ziTzi=1ZiTzi=0\max_z \sum_{i,j} z_i^T \widetilde{K_i}^{-T} K_iK_j^T \widetilde{K_j}^{-1} z_j \\ s.t.\ z_i^T z_i=1 \\ Z_i^Tz_i=0 zmaxi,j∑ziTKi−TKiKjTKj−1zjs.t. ziTzi=1ZiTzi=0

然后可以进一步写为：
max⁡z∑i,jziTPiTKi~−TKiKjTKj~−1Pjzjs.t.ziTzi=1=1\max_z \sum_{i,j} z_i^T P_i^T \widetilde{K_i}^{-T} K_iK_j^T \widetilde{K_j}^{-1} P_jz_j\\ s.t.\ z_i^T z_i=1=1 zmaxi,j∑ziTPiTKi−TKiKjTKj−1Pjzjs.t. ziTzi=1=1

这是因为Pizi=zi−ZiZiTzi=ziP_iz_i=z_i-Z_iZ_i^Tz_i=z_iPizi=zi−ZiZiTzi=zi，这么写把一个约束条件考虑了进去。然后文章在这里再加了一项：
max⁡z∑i,jziTPiTKi~−TKiKjTKj~−1Pjzj+∑iziTzis.t.ziTzi=1=1\max_z \sum_{i,j} z_i^T P_i^T \widetilde{K_i}^{-T} K_iK_j^T \widetilde{K_j}^{-1} P_jz_j+\sum_iz_i^Tz_i \\ s.t.\ z_i^T z_i=1=1 zmaxi,j∑ziTPiTKi−TKiKjTKj−1Pjzj+i∑ziTzis.t. ziTzi=1=1

进一步写为：
max⁡zzTAzs.t.ziTzi=1\max_z z^TAz \\ s.t.\ z_i^T z_i=1 zmaxzTAzs.t. ziTzi=1

个人觉得加的这一项不是很必要，可能是因为核函数的缘故才要这么加，难不成文章认为yiTKiKiTyi=0y_i^TK_iK_i^Ty_i=0yiTKiKiTyi=0吗？不可思议。
AAA是分块矩阵，定义为：
Aij={PiTKi~−TKiKjTKj~−1Pjfor:i≠jIfor:i=jA_{ij}=\left\{ \begin{array}{ccc} P_i^T \widetilde{K_i}^{-T} K_iK_j^T \widetilde{K_j}^{-1} P_j & \ & for: i\neq j \\ I & \ & for: i=j \end{array} \right. Aij={PiTKi−TKiKjTKj−1PjI for:i=jfor:i=j

然后文章里证明了一下AAA矩阵是正定的，这里不写了。因此可以用Horst算法，如下。

真的真的很费劲。理解Horst算法就用d=1d=1d=1就好了。

总结

MCCA在使用上没有什么可用性，也不是这个领域比较好的算法，但是读完这篇论文，我们至少知道了多元特征值问题可以用Horst算法解，并能写出这种名为“MCCA”的算法的代码了。令我有点惊讶的是，直接在网上搜Horst算法没搜出来什么有用的东西，但是一搜多元特征值问题就看到有硕士论文是做这个Horst算法改进的。这个故事告诉我们，在别的领域其实研究的比较多的，已经有经典解的问题，只是因为我们的不熟悉而让自己被卡了壳。所以要多学习，多了解啊。

参考文献

[1] Rupnik J, Shawe-Taylor J. Multi-view canonical correlation analysis[C]//Conference on Data Mining and Data Warehouses (SiKDD 2010). 2010: 1-4.