典型相关分析/典型关联分析 Canonical Correlation Analysis (CCA) 详解过程

我将从以下这几个方面来全面、通俗地为大家介绍典型相关分析 Canonical Correlation Analysis (以下简称，CCA)，及其kernel版本的CCA，简称为KCCA。如果有任何问题，欢迎评论，互相学习讨论！

Why CCA emerges?

为什么CCA会出现呢？因为存在问题：假如我们有两组变量，$X=[x_1,x_2,...,x_m]\in {\mathbb{R}}^{d_m\times m}$（即，$X$有$m$维特征，每维特征有$d_m$个采样点），$Y=[y_1,y_2,...,y_n]\in {\mathbb{R}}^{d_n\times n}$（$Y$有$n$维特征，每个特征有$d_n$个的采样点），如何计算它们的相关性呢？

一般来说，我们首先想到，计算这两组变量的协方差矩阵，$cov(X,Y)$，其中矩阵的每一个元素表示了$X$的第$i（i=1,...,m）$维特征$x_i$与$Y$的第$j（j=1,...,n）$维特征$y_j$之间的相关性。但是，这个做法的缺点是，我们只分析了不同组变量间特征与特征的关系，而忽略了变量内部的关系（也就是，$X$或$Y$内部的特征间的关系），同时也忽略了两组变量的整体相关性（也就是，由$m$维特征构成的$X$与由$n$维特征构成的$Y$之间的关系）。

所以CCA就诞生了：CCA考察两个整体之间的关系，将整体表示成为两组变量各自特征的线性组合。也就是，$Xa$与$Yb$之间的关系。其中向量$a=[a_1,a_2,...,a_m{]}^T$，$b=[b_1,b_2,...,b_m{]}^T$。
$u=Xa=a_1{x}_1+a_2{x}_2+...+a_m{x}_m$
$v=Yb=b_1{y}_1+b_2{y}_2+...+b_n{y}_n$

What is CCA?

CCA是一种将每组变量视为各自特征线性组合形式的整体，通过最大化两个整体的相关系数，从而获得连接系数（$a,b$）、典型相关变量（$u,v$）和整体相关系数$\rho$的一种方法，可以被视为一种降维技术（假设要把两组变量降成一维则只需要选择最大整体相关系数所对应的典型相关变量来表示各自的整体，同理降成二维则增加次大整体相关系数所对应的典型相关变量，以此类推）。目标函数为：

$$\rho =\underset{a,b}{\max }corr(u,v)=\underset{a,b}{\max }\frac{u^{T}v}{\Vert u\Vert \cdot \Vert v\Vert }=\underset{a,b}{\max }\frac{a^T{X}^{T}Yb}{\sqrt{a^T{X}^{T}Xa}\cdot \sqrt{b^T{Y}^{T}Yb}}=\underset{a,b}{\max }\frac{a^T{C}_{\mathbf{x}\mathbf{y}}b}{\sqrt{a^T{C}_{\mathbf{x}\mathbf{x}}a}\cdot \sqrt{b^T{C}_{\mathbf{y}\mathbf{y}}b}}\text{\hspace{1em}(1)}$$
其中，$C_{\mathbf{x}\mathbf{y}}$代表协方差矩阵或方差矩阵。需要说明的是，$X$和$Y$内部的特征假设已经进行了去均值处理，所以相关系数：${\rho }_{x,y}=\frac{(x-\overline{x}{)}^T(y-\overline{y})}{\Vert x-\overline{x}\Vert \cdot \Vert y-\overline{y}\Vert }=\frac{x^{T}y}{\Vert x\Vert \cdot \Vert y\Vert }$，其中$\bar{x}$表示$X$整体的特征均值。

How CCA works?

可以参考$[1]$中的 2.CCA表示及求解。需要提出来的是，与线性判别分析（LDA）类似的是，为了最大化相关系数，在求解时，固定分母而最大化分子。所以，CCA的目标函数可转为优化问题：

• Maximize $a^T{X}^{T}Yb$
• Subject to: $a^T{X}^{T}Xa=1$, $b^T{Y}^{T}Yb=1$

CCA Example

可以参考$[1]$中的 3.CCA计算例子，里面说地通俗易懂。
（原来编辑好的文字，一不小心刷新页面未保存，结果都没了…所以让我们把重点放在kernel版本的CCA上吧。）

Why KCCA emerges?

KCCA的出现原因，大家都能想到：就是因为当我们的数据（两组变量）之间关系是非线性的时候，直接用CCA去寻找每组变量的线性组合，效果是不好的。扩展到核空间下，具有generality的性质。（以下内容主要参考$[2]$，如有需要，可阅读原文。）

What is KCCA？

KCCA是强大的统计工具，能特征化两组多维变量的非线性关系。它首先将数据映射到一个更高维的特征空间中，
$$\varphi :\mathbf{x}=({\mathbf{x}}_1,...,{\mathbf{x}}_n)⟼\varphi (\mathbf{x})=({\varphi }_1(\mathbf{x}),...,{\varphi }_N(\mathbf{x}))(n<N)$$
这里$\mathbf{x}$表示一个具有$n$维的变量集，在经过映射函数$\varphi$之后，该变量集由$n$维变成$N$维（更高维）。其实我们可以理解$N$维为无穷维，所以$\varphi$函数并不可求。为了方便理解和避免混淆，下文将用上标的形式来表示样本，i.e.，${\mathbf{x}}^{(1)}=({\mathbf{x}}_1^{(1)},{\mathbf{x}}_2^{(1)},...,{\mathbf{x}}_n^{(1)})$表示原空间下的第1个$n$维的样本。

根据协方差矩阵的定义，假设$X,Y$表示两个具有$m$个样本，$N$个维度的数据集（$m\times N$），我们得到，
$$C_{\mathbf{x}\mathbf{x}}=X^{T}X,C_{\mathbf{x}\mathbf{y}}=X^{T}Y$$

根据CCA的目标函数，我们依旧可以写目标函数（最大化相关系数）为：
$$\rho =\underset{a,b}{\max }\frac{a^T{C}_{\mathbf{x}\mathbf{y}}b}{\sqrt{a^T{C}_{\mathbf{x}\mathbf{x}}a}\sqrt{b^T{C}_{\mathbf{y}\mathbf{y}}b}}=\underset{a,b}{\max }\frac{a^T{X}^{T}Yb}{\sqrt{a^T{X}^{T}Xa}\sqrt{b^T{Y}^{T}Yb}}\text{\hspace{1em}(1)}$$

$a^T$表示转置。但是我们这里发现，$C_{\mathbf{x}\mathbf{y}}$等是$N\times N$的，$a,b$是$N\times 1$的，所以在核空间里不可求。为此，有一种转换的方法可以解决，就是将我们需要求解的$N$维$a,b$转成一个样本量维度的向量，即，$m\times 1$。这种思想和kernel脊回归里求解连接参数$w$的思想一致，只是kernel脊回归里借助的是一种不等式公式，而KCCA这里，是直接的变换：

$$a=X^T\alpha ,b=Y^T\beta$$
这时，$\alpha ,\beta$是$m\times 1$的。所以目标函数转换为，
$$\rho =\underset{\alpha ,\beta }{\max }\frac{{\alpha }^{T}X{X}^{T}Y{Y}^T\beta }{\sqrt{{\alpha }^{T}X{X}^{T}X{X}^T\alpha \cdot {\beta }^{T}Y{Y}^{T}Y{Y}^T\beta }}.\text{\hspace{1em}(2)}$$

令$K_x=X{X}^T,K_y=Y{Y}^T$，即，

$$\rho =\underset{\alpha ,\beta }{\max }\frac{{\alpha }^T{K}_x{K}_y\beta }{\sqrt{{\alpha }^T{K}_x^2\alpha \cdot {\beta }^T{K}_y^2\beta }}.\text{\hspace{1em}(3)}$$

现在我们发现，该目标函数写成了对偶的形式。观察公式(1)和(3)，我们可以发现，$\alpha ,\beta$的rescale并没有影响到目标函数的形式。

所以，KCCA的优化问题与CCA相同，都是希望找到两组最优系数$\alpha ,\beta$，使得两组变量的相关性最大。

• maximize ${\alpha }^T{K}_x{K}_y\beta$,
• subject to: ${\alpha }^T{K}_x^2\alpha =1,{\beta }^T{K}_y^2\beta =1$。

How KCCA works?

根据目标函数(3)以及它的约束条件，对应的Lagrangian为：
$$L(\lambda ,\alpha ,\beta )={\alpha }^T{K}_x{K}_y\beta -\frac{{\lambda }_{\alpha }}{2}({\alpha }^T{K}_x^2\alpha -1)-\frac{{\lambda }_{\beta }}{2}({\beta }^T{K}_y^2\beta -1),$$
分别对$\alpha ,\beta$求偏导，得

$$\frac{\partial L}{\partial \alpha }=K_x{K}_y\beta -{\lambda }_{\alpha }{K}_x^2\alpha =0\text{\hspace{1em}(4)}$$
$$\frac{\partial L}{\partial \beta }=K_y{K}_x\alpha -{\lambda }_{\beta }{K}_y^2\beta =0\text{\hspace{1em}(5)}$$

公式(4)左乘${\alpha }^T$，公式(5)左乘${\beta }^T$，再两式相减得到，
$$\begin{array}{rl}0=& {\alpha }^T{K}_x{K}_y\beta -{\alpha }^T{\lambda }_{\alpha }{K}_x^2\alpha -{\beta }^T{K}_y{K}_x\alpha +{\beta }^T{\lambda }_{\beta }{K}_y^2\beta \\ =& {\lambda }_{\beta }{\beta }^T{K}_y^2\beta -{\lambda }_{\alpha }{\alpha }^T{K}_x^2\alpha \end{array}$$

由优化目标的条件可知，核矩阵和系数不可能为0，所以${\lambda }_{\beta }-{\lambda }_{\alpha }=0$。让$\lambda ={\lambda }_{\beta }={\lambda }_{\alpha }$，结合公式(4)和(5)，得（假设$K_x,K_y$都是可逆的），
$$\begin{array}{rl}\beta & =\frac{K_y^{-1}{K}_y^{-1}{K}_y{K}_x\alpha }{\lambda }\\ & =\frac{K_y^{-1}{K}_x\alpha }{\lambda }\end{array}$$

$$K_x{K}_y{K}_y^{-1}{K}_x\alpha -{\lambda }^2{K}_x^2\alpha =0$$

$$K_x{K}_x\alpha -{\lambda }^2{K}_x^2\alpha =0$$

$$I\alpha ={\lambda }^2\alpha \text{\hspace{1em}(6)}$$

这是个求解特征根、特征向量的问题。我们乍一看，发现当$\lambda =1$的时候，等式(6)恒成立。所以对于$\alpha$来说，我们可以自由取，即，可以取$\alpha$等于一个$m$维的单位向量，此时，$\beta =\frac{1}{\lambda }{K}_y^{-1}{K}_x$。所以，总的来说，只要$K_x,K_y$是可逆的，我们可以得到完美的系数。但是，现实往往不是这么美好的，也就是这两个核矩阵不可逆的情况常常存在。所以，单纯地用KCCA去寻找合适的$\alpha ,\beta$使得相关性最大，是得不到好结果的。那么接下来，我们将讲解这类不可逆的问题是如何得到解决的。

Computational issues?

这部分除了介绍能解决不可逆问题的Regularization方法外（based on Partial Least Squares），也介绍解决样本量过大导致计算量过大的问题的方法（Gram-Schmidt orthogonolisation/incomplete Cholesky decomposition）。

Regularization

（1）在相关系数的分母的每一项方差项内，加入正则项，如$k||\alpha |{|}^2$，
$$\begin{gathered} \rho =\underset{\alpha ,\beta }{\max }\frac{{\alpha }^T{K}_x{K}_y\beta }{\sqrt{({\alpha }^T{K}_x^2\alpha +k||\alpha |{|}^2)\dot{(}{\beta }^T{K}_y^2\beta +k||\beta |{|}^2)}} \\ =\underset{\alpha ,\beta }{\max }\frac{{\alpha }^T{K}_x{K}_y\beta }{\sqrt{({\alpha }^T{K}_x^2\alpha +k{\alpha }^T{K}_x\alpha )\dot{(}{\beta }^T{K}_y^2\beta +k{\beta }^T{K}_y\beta )}}. \end{gathered}$$

以上第二个等号成立的原因是：我们发现将正则项调整一下，并不会改变本身的优化问题。也就是，调整正则项前后，我们对Lagrangain等式进行求偏导，偏导函数是不变的。所以我们的优化问题就变成了如上所示，且需要服从条件：
subject to: $({\alpha }^T{K}_x^2\alpha +k{\alpha }^T{K}_x\alpha )=1,({\beta }^T{K}_y^2\beta +k{\beta }^T{K}_y\beta )=1$.
所以，根据前面求解的过程，我们可以把最后结果化为以下求解特征问题：
$$(K_x+kI{)}^{-1}{K}_y(K_y+kI{)}^{-1}{K}_x\alpha ={\lambda }^2\alpha$$.
具体的求解过程类似，可参见文献$[2]$。

（2）在KCCA的Lagranrian等式后加上二次正则化项，即，$$\begin{gathered} L_{new}=L+\frac{\eta }{2}(||\alpha |{|}^2+||{\beta }^2||) \\ ={\alpha }^T{K}_x{K}_y\beta -\frac{{\lambda }_{\alpha }}{2}({\alpha }^T{K}_x^2\alpha -1)-\frac{{\lambda }_{\beta }}{2}({\beta }^T{K}_y^2\beta -1)+\frac{\eta }{2}(||\alpha |{|}^2+||{\beta }^2||). \end{gathered}$$
根据同样的方法，该优化为一定可解。

Cholesky Decomposition

当训练集很大时，该问题的计算量会变得很大，而且准确率也会受到影响。为了解决这个问题，我们可以使用Partial Gram-Schmidt orthogonolisation的方法（等价于incomplete Cholesky Decomposition）来降低核矩阵的维度。以下为介绍（主要参考文献[2]）：

• 首先，我们要先知道矩阵分解中的三角分解方法：它可以将一个矩阵分解为一个单位下三角 和 一个上三角矩阵的乘积$A=LU$，且有时候可以是和转置矩阵的乘积$A=L{L}^T$。那什么时候可以是转置矩阵呢？下面的定理告诉你：
  定理1：令A是对称的。如果分解$A=LU$是可能的，那么$L{L}^T=A$。
  定理2：令A是对称且正定的。那么，$A$可以被分解成：$L{L}^T=A$。
  三角分解的用途有很多，如解线性方程、求反矩阵或计算行列式等。如果还想深入了解它，可以自行在网上搜索。

三角分解是一种完整的矩阵分解方法。然后，一个核矩阵的完整分解（complete Cholesky Decomposition）的计算量也不小，而且难以处理现实生活中的数据。所以，以下引入Incomplete Cholesky Decomposition。

• Incomplete Cholesky Decomposition（ICD）详细算法可以参见文献[2]。主要思想是：每一次选择核矩阵的一列，通过贪婪最大化估计误差的下界来选择一列。Partial Gram-Schmidt orthogonolisation（PGSO）的方法也可以从文献[2]中得到参考，ICD与PGSO本质上是等价的，ICD是PGSO的对偶实现。

所以基于以上的思想，我们可以有基于 PGSO的KCCA、基于PGSO的KCCA Regularisation.

Q&A

1. 如果有多组变量/集合$X,Y,Z...$，怎么办？如何衡量多组变量/集合的关系？
   A: 有一种方法称为“Generalization of the Canonical Correlation”，参考$[2]$。该方法使得两两集合的距离之和最小。

2. 如果协方差矩阵不可逆时，怎么办？
   A: Regularization的方法，具体有两类方法：（1）在相关系数的分母的每一项方差项内，加入正则项；（2）在KCCA的Lagranrian等式后加上二次正则化项。具体见上文的Regularization部分。

3. 如果数据样本量太大，导致算法效率低，怎么办？
   A: 使用Cholesky decomposition、Incomplete Cholesky decomposition、Partial Gram-Schmidt Orthogonolisation等

4. 可能有人会问：如何用CCA或KCCA来做预测？
   A: 我们可以这样思考：利用CCA或KCCA，可以得到典型相关系数$a,b,\alpha ,\beta$。那么在给定一个新的样本$x_{n}ew$之后，可以利用得到的典型相关系数计算出$u_{n}ew$；再从所有的$u$中寻找离$u_{new}$最近的几个$u^{\ast }$对应的$x^{\ast }$；再得到$x^{\ast }$对应的$y^{\ast }$，计算$y^{\ast }$每一维度的均值，作为新样本$x_{new}$的预测值。这个做法参考博客[1]的理解，具体实验效果未验证过。

$[1]$ http://www.cnblogs.com/jerrylead/archive/2011/06/20/2085491.html
$[2]$ Hardoon D R, Szedmak S, Shawe-Taylor J. Canonical correlation analysis: An overview with application to learning methods[J]. Neural computation, 2004, 16(12): 2639-2664.