多视角子空间学习系列之 CCA 典型相关分析
多视角学习（Multi-view learning）是陶大成提出的一个研究方向。我们都知道，在机器学习中样本可以用不同的特征（如图片可以用LBP、SIFT等特征）来表示，或者从不同的角度进行观察（如从前后左右观察一个对象），甚至是采用不同的传感器来观测（如RGB摄像头，Depth摄像头），这些不同的观测手段、角度或特征在多视角学习中称为“视角”。多视角学习通过对不同的视角进行统一分析和研究，希望能够得到更好的分类或聚类等等的效果。
陶大成把现有的可用于多视角学习的方法分为三类：（1）Co-training 联合训练 （2）Multiple Kernel Learning 多核学习 （3）Subspace Learning 子空间学习。 我所研究的正是Subspace Learning这个方向。
每个样本可以看做是高维空间上分布的一个点，每个视角所有的样本的分布构成一个样本空间，子空间学习（Subspace Learning）认为这些样本空间存在一个潜藏的公共子空间，各视角的各样本在这个公共子空间中都有一个投影，或者叫做表示。子空间学习的目标就是寻找到这个公共子空间，并让各样本在其中的表示具有更好的某些性质，或者保持原始分布的某些性质。如果子空间的维度低于原始样本空间的维度，就产生了降维的问题，因此子空间学习与多视角降维几乎是同一个问题。
子空间学习要面对的主要问题是，不同的视角的数据维度可能是不相同的，这使得跨视角的度量变得困难。
从CCA开始我写博客介绍一些子空间学习方法，我在这个领域的研究是从CCA开始的，我对CCA的学习也是因为当时要研究这个领域。
优化目标CCA（Canonical Correlation Analysis）典型相关分析是最经典的子空间学习方法。它是一种线性无监督的方法，能处理视角数量是两个。给定两个做过中心化的视角X1∈RD1×n,X2∈RD2×nX_1\in \mathbb{R}^{D_1\times n},X_2\in \mathbb{R}^{D_2\times n}X1∈RD1×n,X2∈RD2×n的样本，其中D1,D2D_1,D_2D1,D2表示二者的维度，两者未必相等，两个视角的样本数量均为nnn，一列为一个样本。CCA希望找到两个线性的转换矩阵W1∈RD1×d,W2∈RD2×d,d≤min(D1,D2)W_1\in \mathbb{R}^{D_1\times d},W_2\in \mathbb{R}^{D_2\times d},d\leq min(D_1,D_2)W1∈RD1×d,W2∈RD2×d,d≤min(D1,D2)使得W1TX1W_1^TX_1W1TX1与W2TX2W_2^TX_2W2TX2的相关系数最大。这就是CCA的优化目标和求解目标。
来复习一下相关系数的公式，如果a,ba,ba,b为两列随机变量，a,ba,ba,b的相关系数r(a,b)r(a,b)r(a,b)定义为
r(a,b)=Cov(a,b)Var(a)Var(b)Cov(a,b)=E[a−E(a)]E[b−E(b)]=E(ab)−E(a)E(b)
r(a,b)=\frac{Cov(a,b)}{\sqrt{Var(a)Var(b)}} \\
Cov(a,b)=E[a-E(a)]E[b-E(b)]=E(ab)-E(a)E(b)
r(a,b)=Var(a)Var(b)Cov(a,b)Cov(a,b)=E[a−E(a)]E[b−E(b)]=E(ab)−E(a)E(b)
CovCovCov是协方差，VarVarVar是方差。若Y1=W1TX1,Y2=W2TX2Y_1=W_1^TX_1,Y_2=W_2^TX_2Y1=W1TX1,Y2=W2TX2，CCA的目标就是使二者的相关系数最大化，Y1∈Rd×n,Y2∈Rd×nY_1\in \mathbb{R}^{d\times n},Y_2\in \mathbb{R}^{d\times n}Y1∈Rd×n,Y2∈Rd×n就是X1,X2X_1,X_2X1,X2在公共子空间的表示。
d=1d=1d=1的Lagrangian乘子法解若d=1d=1d=1，W1,W2W_1,W_2W1,W2都为向量，我们改用w1,w2w_1,w_2w1,w2表示。此时：
Cov(w1TX1,w2TX2)=E(w1TX1X2Tw2)−E(w1TX1)E(w2TX2)
Cov(w_1^TX_1,w_2^TX_2)=E(w_1^TX_1X_2^Tw_2)-E(w_1^TX_1)E(w_2^TX_2)
Cov(w1TX1,w2TX2)=E(w1TX1X2Tw2)−E(w1TX1)E(w2TX2)
这里E(w1TX1X2Tw2)E(w_1^TX_1X_2^Tw_2)E(w1TX1X2Tw2)是一个数字，其期望就是本身；E(w1TX1),E(w2TX2)E(w_1^TX_1),E(w_2^TX_2)E(w1TX1),E(w2TX2)是向量，其期望是平均值，而我们认为X1,X2X_1,X_2X1,X2都已经做过中心化了，即X1,X2X_1,X_2X1,X2的行均值均为0，因此E(w1TX1),E(w2TX2)E(w_1^TX_1),E(w_2^TX_2)E(w1TX1),E(w2TX2)都为000。
方差是咋算的呢，以w1TX1w_1^TX_1w1TX1为例，同样是因为已经做过中心化了：
Var(w1TX1)=1n−1w1TX1X1Tw1
Var(w_1^TX_1)=\frac{1}{n-1}w_1^TX_1X_1^Tw_1
Var(w1TX1)=n−11w1TX1X1Tw1
此时优化目标变成了：
maxw1,w2w1TX1X2Tw2w1TX1X1Tw1n−1w2TX2X2Tw2n−1
\max_{w_1,w_2} \frac{w_1^TX_1X_2^Tw_2}{ \sqrt{\frac{w_1^TX_1X_1^Tw_1}{n-1} \frac{w_2^TX_2X_2^Tw_2}{n-1} }}
w1,w2maxn−1w1TX1X1Tw1n−1w2TX2X2Tw2w1TX1X2Tw2
注意到同一对w1,w2w_1,w_2w1,w2，如果同时扩大ccc倍变成cw1,cw2cw_1,cw_2cw1,cw2，上式的值不变，因此为了求解方便，也为了上式有唯一解，CCA添加限制条件w1TX1X1Tw1=1,w2TX2X2Tw2=1w_1^TX_1X_1^Tw_1=1,w_2^TX_2X_2^Tw_2=1w1TX1X1Tw1=1,w2TX2X2Tw2=1。这样总的求解目标就等价于：
maxw1,w2w1TX1X2Tw2s.t. w1TX1X1Tw1=1,w2TX2X2Tw2=1
\max_{w_1,w_2}w_1^TX_1X_2^Tw_2 \\
s.t.\ w_1^TX_1X_1^Tw_1=1,w_2^TX_2X_2^Tw_2=1
w1,w2maxw1TX1X2Tw2s.t. w1TX1X1Tw1=1,w2TX2X2Tw2=1
一般论文里列出来的式子就是这个式子。虽然曲折，但是每一步都有理有据。这是一个有约束的凸优化问题（为什么是凸函数呢？因为二阶导是0矩阵，是半正定的），可以用Lagrangian乘子法来解：
L(w1,w2,λ1,λ2)=w1TX1X2Tw2+λ1(1−w1TX1X1Tw1)+λ2(1−w2TX2X2Tw2)
L(w_1,w_2,\lambda_1,\lambda_2)=w_1^TX_1X_2^Tw_2+\lambda_1(1-w_1^TX_1X_1^Tw_1)+\lambda_2(1-w_2^TX_2X_2^Tw_2)
L(w1,w2,λ1,λ2)=w1TX1X2Tw2+λ1(1−w1TX1X1Tw1)+λ2(1−w2TX2X2Tw2)
对四个变量求偏导并令结果为0：
∂∂w1L(w1,w2,λ1,λ2)=X1X2Tw2−2λ1X1X1Tw1=0X1X2Tw2=2λ1X1X1Tw1∂∂w2L(w1,w2,λ1,λ2)=X2X1Tw1−2λ2X2X2Tw2=0X2X1Tw1=2λ2X2X2Tw2∂∂λ1L(w1,w2,λ1,λ2)=1−w1TX1X1Tw1=0w1TX1X1Tw1=1∂∂λ2L(w1,w2,λ1,λ2)=1−w2TX2X2Tw2=0w2TX2X2Tw2=1
\frac{\partial }{\partial w_1}L(w_1,w_2,\lambda_1,\lambda_2)=X_1X_2^Tw_2-2\lambda_1 X_1X_1^Tw_1=0 \\
X_1X_2^Tw_2=2\lambda_1 X_1X_1^Tw_1 \\
\frac{\partial }{\partial w_2}L(w_1,w_2,\lambda_1,\lambda_2)=X_2X_1^Tw_1-2\lambda_2 X_2X_2^Tw_2=0 \\
X_2X_1^Tw_1=2\lambda_2 X_2X_2^Tw_2 \\
\frac{\partial }{\partial \lambda_1}L(w_1,w_2,\lambda_1,\lambda_2) =1-w_1^TX_1X_1^Tw_1=0\\
w_1^TX_1X_1^Tw_1=1 \\
\frac{\partial }{\partial \lambda_2}L(w_1,w_2,\lambda_1,\lambda_2)=1-w_2^TX_2X_2^Tw_2=0 \\
w_2^TX_2X_2^Tw_2=1
∂w1∂L(w1,w2,λ1,λ2)=X1X2Tw2−2λ1X1X1Tw1=0X1X2Tw2=2λ1X1X1Tw1∂w2∂L(w1,w2,λ1,λ2)=X2X1Tw1−2λ2X2X2Tw2=0X2X1Tw1=2λ2X2X2Tw2∂λ1∂L(w1,w2,λ1,λ2)=1−w1TX1X1Tw1=0w1TX1X1Tw1=1∂λ2∂L(w1,w2,λ1,λ2)=1−w2TX2X2Tw2=0w2TX2X2Tw2=1
有用的主要就这两个，剩下俩是约束条件：
X1X2Tw2=2λ1X1X1Tw1X2X1Tw1=2λ2X2X2Tw2
X_1X_2^Tw_2=2\lambda_1 X_1X_1^Tw_1 \\
X_2X_1^Tw_1=2\lambda_2 X_2X_2^Tw_2
X1X2Tw2=2λ1X1X1Tw1X2X1Tw1=2λ2X2X2Tw2
做如下推导：
w1TX1X2Tw2=2λ1w1TX1X1Tw1=2λ1w2TX2X1Tw1=2λ2w2TX2X2Tw2=2λ2∴λ1T=λ2,λ1=λ2
w_1^TX_1X_2^Tw_2=2\lambda_1 w_1^TX_1X_1^Tw_1=2\lambda_1 \\
w_2^T X_2X_1^Tw_1=2\lambda_2 w_2^T X_2X_2^Tw_2 =2\lambda_2 \\
\therefore \lambda_1^T=\lambda_2,\lambda_1=\lambda_2
w1TX1X2Tw2=2λ1w1TX1X1Tw1=2λ1w2TX2X1Tw1=2λ2w2TX2X2Tw2=2λ2∴λ1T=λ2,λ1=λ2
令λ=λ1=λ2\lambda=\lambda_1=\lambda_2λ=λ1=λ2：
(X1X1T)−1X1X2Tw2=2λw1(X2X2T)−1X2X1Tw1=2λw2
(X_1X_1^T)^{-1}X_1X_2^Tw_2=2\lambda w_1 \\
(X_2X_2^T)^{-1}X_2X_1^Tw_1=2\lambda w_2
(X1X1T)−1X1X2Tw2=2λw1(X2X2T)−1X2X1Tw1=2λw2
进一步得到：
(X1X1T)−1X1X2T(X2X2T)−1X2X1Tw1=4λ2w1(X2X2T)−1X2X1T(X1X1T)−1X1X2Tw2=4λ2w2
(X_1X_1^T)^{-1}X_1X_2^T(X_2X_2^T)^{-1}X_2X_1^Tw_1=4\lambda^2 w_1 \\
(X_2X_2^T)^{-1}X_2X_1^T(X_1X_1^T)^{-1}X_1X_2^Tw_2=4\lambda^2 w_2
(X1X1T)−1X1X2T(X2X2T)−1X2X1Tw1=4λ2w1(X2X2T)−1X2X1T(X1X1T)−1X1X2Tw2=4λ2w2
这样就解出来w1,w2w_1,w_2w1,w2了，分别是那两个很复杂的矩阵的特征向量，λ\lambdaλ是特征值。因为是个必要条件，而且前面已经得到w1TX1X2Tw2=2λw_1^TX_1X_2^Tw_2=2\lambdaw1TX1X2Tw2=2λ，因此只要选择最大的特征值对应的特征向量就能使得原式最大化。
d>1d>1d>1的Lagrangian乘子法解若d>1d>1d>1，问题变为：
maxW1,W2tr(W1TX1X2TW2)s.t. W1TX1X1TW1=I,W2TX2X2TW2=1
\max_{W_1,W_2} tr ( W_1^TX_1X_2^TW_2) \\
s.t.\ W_1^TX_1X_1^TW_1=I,W_2^TX_2X_2^TW_2=1
W1,W2maxtr(W1TX1X2TW2)s.t. W1TX1X1TW1=I,W2TX2X2TW2=1
即在各维度上的相关系数之和最大。列Lagrangian乘子法：
L(W1,W2,λ1,λ2)=tr(W1TX1X2TW2)+tr[λ1(I−W1TX1X1TW1)]+tr[λ2(I−W2TX2X2TW2)]
L(W_1,W_2,\lambda_1,\lambda_2)= tr(W_1^TX_1X_2^TW_2)+tr[\lambda_1(I-W_1^TX_1X_1^TW_1)]+tr[\lambda_2(I-W_2^TX_2X_2^TW_2)]
L(W1,W2,λ1,λ2)=tr(W1TX1X2TW2)+tr[λ1(I−W1TX1X1TW1)]+tr[λ2(I−W2TX2X2TW2)]
其中λ1,λ2\lambda_1,\lambda_2λ1,λ2都是对角矩阵。求导并令导数为0：
∂∂W1L(W1,W2,λ1,λ2)=X1X2TW2−2X1X1TW1λ1=0X1X2TW2=2X1X1TW1λ1∂∂W2L(W1,W2,λ1,λ2)=X2X1TW1−2X2X2TW2λ2=0X2X1TW1=2X2X2TW2λ2
\frac{\partial }{\partial W_1}L(W_1,W_2,\lambda_1,\lambda_2)=X_1X_2^TW_2-2X_1X_1^TW_1\lambda_1=0 \\
X_1X_2^TW_2=2X_1X_1^TW_1\lambda_1 \\
\frac{\partial }{\partial W_2}L(W_1,W_2,\lambda_1,\lambda_2) = X_2X_1^TW_1-2X_2X_2^TW_2\lambda_2 =0\\
X_2X_1^TW_1=2X_2X_2^TW_2\lambda_2 \\
∂W1∂L(W1,W2,λ1,λ2)=X1X2TW2−2X1X1TW1λ1=0X1X2TW2=2X1X1TW1λ1∂W2∂L(W1,W2,λ1,λ2)=X2X1TW1−2X2X2TW2λ2=0X2X1TW1=2X2X2TW2λ2
另外两个对λ1,λ2\lambda_1,\lambda_2λ1,λ2求导的我们就不写了，没必要。然后做转化，证明λ1=λ2\lambda_1=\lambda_2λ1=λ2：
W1TX1X2TW2=2W1TX1X1TW1λ1=2λ1W2TX2X1TW1=2W2TX2X2TW2λ2=2λ2∴λ1T=λ2,λ1=λ2
W_1^TX_1X_2^TW_2=2 W_1^TX_1X_1^TW_1\lambda_1=2\lambda_1 \\
W_2^TX_2X_1^TW_1=2W_2^TX_2X_2^TW_2\lambda_2 =2\lambda_2 \\
\therefore \lambda_1^T=\lambda_2,\lambda_1=\lambda_2
W1TX1X2TW2=2W1TX1X1TW1λ1=2λ1W2TX2X1TW1=2W2TX2X2TW2λ2=2λ2∴λ1T=λ2,λ1=λ2
令λ=λ1=λ2\lambda=\lambda_1=\lambda_2λ=λ1=λ2，然后重写一下：
X1X2TW2=2X1X1TW1λX2X1TW1=2X2X2TW2λ
X_1X_2^TW_2=2X_1X_1^TW_1\lambda\\
X_2X_1^TW_1=2X_2X_2^TW_2\lambda
X1X2TW2=2X1X1TW1λX2X1TW1=2X2X2TW2λ
计算：
W1=12(X1X1T)−1(X1X2T)W2λ−1X2X1T(X1X1T)−1X1X2TW2=4(X2X2T)W2λ2(X2X2T)−1X2X1T(X1X1T)−1X1X2TW2=4W2λ2
W_1=\frac{1}{2}(X_1X_1^T)^{-1}(X_1X_2^T)W_2\lambda^{-1} \\
X_2X_1^T(X_1X_1^T)^{-1}X_1X_2^TW_2=4(X_2X_2^T)W_2\lambda^2 \\
(X_2X_2^T)^{-1}X_2X_1^T(X_1X_1^T)^{-1}X_1X_2^TW_2=4W_2\lambda^2
W1=21(X1X1T)−1(X1X2T)W2λ−1X2X1T(X1X1T)−1X1X2TW2=4(X2X2T)W2λ2(X2X2T)−1X2X1T(X1X1T)−1X1X2TW2=4W2λ2
同理：
(X1X1T)−1X1X2T(X2X2T)−1X2X1TW1=4W2λ2
(X_1X_1^T)^{-1}X_1X_2^T(X_2X_2^T)^{-1}X_2X_1^TW_1=4W_2 \lambda^2
(X1X1T)−1X1X2T(X2X2T)−1X2X1TW1=4W2λ2
这就告诉我们W1,W2W_1,W_2W1,W2的每一列都分别是由上面复杂的那两个式子的特征向量构成的，对应的特征值分布在λ2\lambda^2λ2的对角线上。因为这是个必要条件，现在把W1TX1X2TW2=2λW_1^TX_1X_2^TW_2=2\lambdaW1TX1X2TW2=2λ代回求解目标：
tr(W1TX1X2TW2)=tr(2λ)
tr ( W_1^TX_1X_2^TW_2)=tr(2\lambda)
tr(W1TX1X2TW2)=tr(2λ)
为了使求解目标最大化，就要使tr(2λ)tr(2\lambda)tr(2λ)最大，因此我们要选上面那两个复杂表示的最大的特征值对应的ddd个特征向量构成W1,W2W_1,W_2W1,W2。
总结CCA是非常经典的方法，在数据分析、金融等领域应用广泛，CCA之于子空间学习相当于PCA之于降维。
作者：魏之燕
