多视角子空间学习系列之 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−1w1TX1X1Tw1Var(w_1^TX_1)=\frac{1}{n-1}w_1^TX_1X_1^Tw_1 Var(w1TX1)=n−11w1TX1X1Tw1

此时优化目标变成了：
max⁡w1,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。这样总的求解目标就等价于：
max⁡w1,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λ2X2X2Tw2X_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=λ2w_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，问题变为：
max⁡W1,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=λ2W_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λ2W_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之于降维。