典型性相关分析Python代码 典型相关性分析的特点

典型相关分析

 （一）引入

     典型相关分析（Canonical Correlation Analysis）是研究两组变量之间相关关系的一种多元统计方法。他能够揭示出两组变量之间的内在联系。

    我们知道，在一元统计分析中，用相关系数来衡量两个随机变量的线性相关关系，用复相关系数研究一个随机变量与多个随机变量的线性相关关系。然而，这些方法均无法用于研究两组变量之间的相关关系，于是提出了CCA。其基本思想和主成分分析非常相似。首先，在每组变量中寻找出变量的线性组合，使得两组的线性组合之间具有最大的相关系数；然后选取和已经挑选出的这对线性组合不相关的另一对线性组合，并使其相关系数最大，如此下去，直到两组变量的相关性被提取完毕为止。被选出的线性组合配对称为典型变量，它们的相关系数称为典型相关系数。

（二）分析

    设有两组随机变量$\mathbf{X}=(x_1,x_2,\cdots,x_p)^\prime$和$\mathbf{Y}=(y_1,y_2,\cdots,y_q)^\prime$，不妨设$p\leq q$。设第一组变量均值为$\mathbb{E}\mathbf{X}=\mu_1$，方差为$\mathop{Var}(\mathbf{X})=\mathop{cov}(\mathbf{X},\mathbf{X})=\Sigma_{11}$。第二组变量均值为$\mathbb{E}\mathbf{Y}=\mu_2$，方差为$\mathop{Var}(\mathbf{Y})=\mathop{cov}(\mathbf{Y},\mathbf{Y})=\Sigma_{22}$。第一组与第二组变量的协方差矩阵为$\mathop{cov}(\mathbf{X},\mathbf{Y})=\Sigma_{12}=\Sigma_{21}^\prime$。

    分别对两组变量做线性组合：

\begin{align}U &= a_1x_1+a_2x_2+\cdots+a_px_p=a^\prime\mathbf{X}\label{equ:U}\\ V &= b_1x_1+b_2x_2+\cdots+b_qy_q=b^\prime\mathbf{Y}\label{equ:V}\end{align}

所以$U,V$的方差，协方差，相关系数为：

\begin{align}&\mathop{Var}(U)=a^\prime\mathop{cov}(\mathbf{X},\mathbf{X})a=a^\prime\Sigma_{11}a\label{equ:UVar}\\&\mathop{Var}(V)=b^\prime\mathop{cov}(\mathbf{Y},\mathbf{Y})b=b^\prime\Sigma_{22}b\label{VVar}\\&\mathop{cov}(U,V)=a^\prime\mathop{cov}(\mathbf{X},\mathbf{Y})b=a^\prime\Sigma_{12}b\label{equ:UVcov}\\&\rho=\mathop{corr}(U,V)=\frac{a^\prime\Sigma_{12}b}{\sqrt{a^\prime\Sigma_{11}a}\sqrt{b^\prime\Sigma_{12}b}}\label{equ:rho}\end{align}

其中$U,V$称为典型变量，它们之间的相关系数$\rho$称为典型相关系数。

    CCA要解决的问题是，在所有线性组合$U$和$V$中选取典型相关系数最大的那对，即选取$a^{(1)},b^{(1)}$使$U_1=(a^{(1)})^\prime\mathbf{X}$与$V_1=(b^{(1)})^\prime\mathbf{Y}$之间的相关系数最大，这里$(U_1,V_1)$称为第一对典型相关变量；然后在选取$a^{(2)},b^{(2)}$使得$U_1=(a^{(2)})^\prime\mathbf{X},V_2(b^{(2)})^\prime\mathbf{Y}$，在与$U_1,V_1$不相关的情况下，使得$(U_2,V_2)$的相关系数最大，称为第二对典型相关变量；如此继续下去，直到所有分别与$(U_1,V_1),(U_2,V_2),\cdots,(U_{p-1},V_{p-1})$都不相关的线性组合$(U_p,V_p)$为止，此时$p$为$\mathbf{X}$与$\mathbf{Y}$之间的协方差矩阵的秩。由上面的分析可得模型：

\begin{equation}\mathop{\max}\quad\rho=\frac{a^\prime\Sigma_{12}b}{\sqrt{a^\prime\Sigma_{11}a}\sqrt{b^\prime\Sigma_{22}b}}\label{model:original}\end{equation}

由于收缩$U$和$V$的值并不会影响$\rho$，故我们可引入限制条件$a^\prime\Sigma_{11}a=1,b^\prime\Sigma_{22}b=1$将模型转化为：

\begin{align}\mathop{\max}\quad&a^\prime\Sigma_{12}b\nonumber\\\mathop{s.t.}\quad&a^\prime\Sigma_{11}a=1\nonumber\\\quad& b^\prime\Sigma_{22}b=1\label{model:constraint}\end{align}

引入Lagrange乘子：

\begin{equation}L(a,b,\lambda,\nu)=a^\prime\Sigma_{12}b-\frac{\lambda}{2}(a^\prime\Sigma_{11}a-1)-\frac{\nu}{2}(b^\prime\Sigma_{22}b-1)\label{equ:lagrange}\end{equation}

对Lagrange函数\ref{equ:lagrange}求导得：

\begin{equation}\frac{\partial L}{\partial a}=\Sigma_{12}b-\lambda\Sigma_{11}a=0\label{equ:partiala}\end{equation}

\begin{equation}\frac{\partial L}{\partial b}=\Sigma_{21}a-\nu\Sigma_{22}b=0\label{equ:partialb}\end{equation}

将式子\ref{equ:partiala}左乘$a^\prime$，式子\ref{equ:partialb}左乘$b^\prime$得：

\begin{align*}a^\prime\Sigma_{12}b-\lambda a^\prime\Sigma_{11}a &=0\\b^\prime\Sigma_{21}a-\nu b^\prime\Sigma_{22}b &=0\end{align*}

 又因为$(a^\prime\Sigma_{12}b)^\prime=b^\prime\Sigma_{21}a\Longrightarrow \lambda a^\prime\Sigma_{11}a=\nu b^\prime\Sigma_{22}b$。由限制条件知：$\lambda=\nu=\rho=a^\prime\Sigma_{12}b$，即$\lambda$的值就是线性组合$U$和$V$的相关系数。我们重新将式子\ref{equ:partiala}和式子\ref{equ:partialb}写成：

\begin{align}-\lambda\Sigma_{11}a+\Sigma_{12}b&=0\label{equ:partialaNew}\\\Sigma_{21}a-\lambda\Sigma_{22}b&=0\label{equ:partialbNew}\end{align}

然后将式子\ref{equ:partialbNew}左乘$\Sigma_{12}\Sigma_{22}^{-1}$得：

\begin{equation}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}a=\lambda\Sigma_{12}b\label{equ:partialbLeft}\end{equation}

结合式子\ref{equ:partialaNew}得：

\begin{equation}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}a=\lambda^2\Sigma_{11}a\Longrightarrow (\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}-\lambda^2\Sigma_{11})a=0\label{equ:3to4}\end{equation}

同理，将式子\ref{equ:partialaNew}左乘$\Sigma_{21}\Sigma_{11}^{-1}$,并将式子\ref{equ:partialbNew}代入式子\ref{equ:partialaNew}得：

\begin{equation}(\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}-\lambda^2\Sigma_{22})b=0\label{equ:4to3}\end{equation}

将$\Sigma_{11}^{-1}$左乘式子\ref{equ:3to4}，$\Sigma_{22}^{-1}$左乘式子\ref{equ:4to3}得：

\begin{equation}\left\{\begin{array}&(\Sigma_{11}^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}-\lambda^2)a=0\\(\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}-\lambda^2)b=0\end{array}\right.\triangleq\left\{\begin{array}&\mathbf{A}a=\lambda^2a\\\mathbf{B}b=\lambda^2b\end{array}\right.\end{equation}

说明：$\lambda^2$既是矩阵$\mathbf{A}$也是矩阵$\mathbf{B}$的特征值，$a$与$b$分别是对应的特征向量。所以我们的问题转化成求矩阵$\mathbf{A},\mathbf{B}$的最大特征值对应的特征向量，而特征值的平方根$\sqrt{\lambda}$为相关系数，从而求出第一对典型相关变量。

    此时，我们可以得到如下的猜想：是否矩阵$\mathbf{A},\mathbf{B}$的所有非零特征值的平方跟都会是其对应的典型相关系数？接下去，我们来证明如下猜想：

设$\lambda^2_1\geq\lambda_2^2\geq\cdots\geq\lambda_p^2>0$为$\mathbf{A},\mathbf{B}$的特征值（$p$为矩阵$\Sigma_{12}$的秩），其对应的特征向量为$a^{(1)},a^{(2)},\cdots,a^{(p)}$和$b^{(1)},b^{(2)},\cdots,b^{(p)}$，于是得到$p$对线性组合：

\begin{equation}\begin{array}&U_1=(a^{(1)})^\prime\mathbf{X}&U_2=(a^{(2)})^\prime\mathbf{X}&\cdots&U_p=(a^{(p)})^\prime\mathbf{X}\\V_1=(b^{(1)})^\prime\mathbf{Y}&V_2=(b^{(2)})^\prime\mathbf{Y}&\cdots&V_p=(b^{(p)})^\prime\mathbf{Y}\end{array}\end{equation}

可以证明$(U_1,V_1),(U_2,V_2),\cdots,(U_p,V_p)$就是其前$p$对典型变量，$\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_p$为其典型相关系数。

首先，在求出第一对典型变量的基础上求第二对典型变量。由上述分析我们可以知道该模型为：

\begin{align}\mathop{\max}&\quad (a^{(2)})^\prime\Sigma_{12}b^{(2)}\nonumber\\\mathop{s.t.}&\quad (a^{(2)})^\prime\Sigma_{11} a^{(2)}=1\nonumber\\&\quad(b^{(2)})^\prime\Sigma_{22}b^{(2)}=1\nonumber\\&\quad(a^{(1)})^\prime\Sigma_{11}a^{(2)}=0\label{con:1}\\&\quad(b^{(1)})^\prime\Sigma_{22}b^{(2)}=0\label{con:2}\end{align}

其中限制\ref{con:1}和\ref{con:2}是由于第二对典型变量必须与第一对典型变量无关。其Lagrange方程以及相应的偏导为：

\begin{align*}L(a^{(1)},b^{(1)},a^{(2)},b^{(2)}) &=(a^{(2)})^\prime\Sigma_{12}b^{(2)}-\frac{\lambda}{2}((a^{(2)})^\prime\Sigma_{11}a^{(2)}-1)-\frac{\nu}{2}((b^{(2)})^\prime\Sigma_{22}b^{(2)}-1)\\&\quad + \gamma (a^{(1)})^\prime\Sigma_{11}a^{(2)}+\beta(b^{(1)})^\prime\Sigma_{22}b^{(2)}\end{align*}

\begin{align}\frac{\partial L}{\partial a^{(2)}} &=\Sigma_{12}b^{(2)}-\lambda\Sigma_{11}a^{(2)}+\gamma\Sigma_{11}a^{(1)}=0\label{partial:a2}\\\frac{\partial L}{\partial b^{(2)}}&=\Sigma_{21}a^{(2)}-\nu\Sigma_{22}b^{(2)}+\beta\Sigma_{22}b^{(1)}=0\label{partial:b2}\\\frac{\partial L}{\partial a1}&=\gamma\Sigma_{11}a^{(2)}=0\label{partial:a1}\\\frac{\partial L}{\partial b^{(1)}}&=\beta\Sigma_{22}b^{(2)}=0\label{partial:b1}\end{align}

将$(a^{(2)})^\prime$左乘式子\ref{partial:a2}，$(b^{(2)})^\prime$左乘式子\ref{partial:b2}得：

\begin{align}&(a^{(2)})^\prime\Sigma_{12}b^{(2)}-\lambda(a^{(2)})^\prime\Sigma_{11}a^{(2)}+\gamma (a^{(2)})^\prime\Sigma_{11}a^{(1)}=0\label{partial:a2new}\\&(b^{(2)})^\prime\Sigma_{21}a^{(2)}-\nu(b^{(2)})^\prime\Sigma_{22}b^{(2)}+\beta (b^{(2)})^\prime\Sigma_{22}b^{(1)}=0\label{partial:b2new}\end{align}

将$(a^{(1)})^\prime$左乘式子\ref{partial:a1}，$(b^{(2)})^\prime$左乘式子\ref{partial:b1}得：

\begin{align}&(a^{(2)})^\prime\Sigma_{11}a^{(1)}=0\label{partial:a1new}\\&(b^{(2)})^\prime\Sigma_{22}b^{(2)}=0\label{partial:b1new}\end{align}

将式子\ref{partial:a1new}和式子\ref{partial:b1new}代入式子\ref{partial:a2new}和\ref{partial:b2new}得：

\begin{align}&(a^{(2)})^\prime\Sigma_{12}b^{(2)}-\lambda(a^{(2)})^\prime\Sigma_{11}a^{(2)}=0\label{partial:a1a2}\\&(b^{(2)})^\prime\Sigma_{21}a^{(2)}-\nu(b^{(2)})^\prime\Sigma_{22}b^{(2)}=0\label{partial:b1b2}\end{align}

其中式子\ref{partial:a1a2}和式子\ref{partial:b1b2}与式子\ref{equ:partiala}和式子\ref{equ:partialb}有相同的形式，只是$a,b$换成$a^{(2)},b^{(2)}$，故同样可以得到：

\begin{equation}\left\{\begin{array}&\mathbf{A}a^{(2)}=\lambda^2a^{(2)}\\\mathbf{B}b^{(2)}=\lambda^2b^{(2)}\end{array}\right.\end{equation}

此时$a^{(2)}\neq a,b^{(2)}\neq b$，否则不满足限制\ref{con:1}和\ref{con:2}，所以最大值为第二大特征值。以此类推，我们即可证明上述猜想。

注意：我们在求解上述普通特征值方程$\mathbf{A}a=\lambda^2a$时，由于$\mathbf{A}=\Sigma_{11}^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$而求逆过程的计算量大，精度低，故我们可以将其中的对称矩阵$\Sigma_{11}$进行Cholesky分解，即$\Sigma_{11}=\mathbf{R_1}\mathbf{R_1}^\prime$，其中$\mathbf{R_1}$为下三角矩阵。于是方程可化为对称矩阵的求特征值问题：$\mathbf{R_1}^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}(\mathbf{R_1}^{-1})^\prime u_x=\lambda^2 u_x$，其中$u_x=\mathbf{R_1}^\prime a$。