典型相关分析
(一)引入
典型相关分析(Canonical Correlation Analysis)是研究两组变量之间相关关系的一种多元统计方法。他能够揭示出两组变量之间的内在联系。
我们知道,在一元统计分析中,用相关系数来衡量两个随机变量的线性相关关系,用复相关系数研究一个随机变量与多个随机变量的线性相关关系。然而,这些方法均无法用于研究两组变量之间的相关关系,于是提出了CCA。其基本思想和主成分分析非常相似。首先,在每组变量中寻找出变量的线性组合,使得两组的线性组合之间具有最大的相关系数;然后选取和已经挑选出的这对线性组合不相关的另一对线性组合,并使其相关系数最大,如此下去,直到两组变量的相关性被提取完毕为止。被选出的线性组合配对称为典型变量,它们的相关系数称为典型相关系数。
(二)分析
设有两组随机变量$mathbf{X}=(x_1,x_2,cdots,x_p)^prime$和$mathbf{Y}=(y_1,y_2,cdots,y_q)^prime$,不妨设$pleq 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$。
分别对两组变量做线性组合:
egin{align}U &= a_1x_1+a_2x_2+cdots+a_px_p=a^primemathbf{X}label{equ:U}\ V &= b_1x_1+b_2x_2+cdots+b_qy_q=b^primemathbf{Y}label{equ:V}end{align}
所以$U,V$的方差,协方差,相关系数为:
egin{align}&mathop{Var}(U)=a^primemathop{cov}(mathbf{X},mathbf{X})a=a^primeSigma_{11}alabel{equ:UVar}\&mathop{Var}(V)=b^primemathop{cov}(mathbf{Y},mathbf{Y})b=b^primeSigma_{22}blabel{VVar}\&mathop{cov}(U,V)=a^primemathop{cov}(mathbf{X},mathbf{Y})b=a^primeSigma_{12}blabel{equ:UVcov}\& ho=mathop{corr}(U,V)=frac{a^primeSigma_{12}b}{sqrt{a^primeSigma_{11}a}sqrt{b^primeSigma_{12}b}}label{equ:rho}end{align}
其中$U,V$称为典型变量,它们之间的相关系数$ ho$称为典型相关系数。
CCA要解决的问题是,在所有线性组合$U$和$V$中选取典型相关系数最大的那对,即选取$a^{(1)},b^{(1)}$使$U_1=(a^{(1)})^primemathbf{X}$与$V_1=(b^{(1)})^primemathbf{Y}$之间的相关系数最大,这里$(U_1,V_1)$称为第一对典型相关变量;然后在选取$a^{(2)},b^{(2)}$使得$U_1=(a^{(2)})^primemathbf{X},V_2(b^{(2)})^primemathbf{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}$之间的协方差矩阵的秩。由上面的分析可得模型:
egin{equation}mathop{max}quad ho=frac{a^primeSigma_{12}b}{sqrt{a^primeSigma_{11}a}sqrt{b^primeSigma_{22}b}}label{model:original}end{equation}
由于收缩$U$和$V$的值并不会影响$ ho$,故我们可引入限制条件$a^primeSigma_{11}a=1,b^primeSigma_{22}b=1$将模型转化为:
egin{align}mathop{max}quad&a^primeSigma_{12}b onumber\mathop{s.t.}quad&a^primeSigma_{11}a=1 onumber\quad& b^primeSigma_{22}b=1label{model:constraint}end{align}
引入Lagrange乘子:
egin{equation}L(a,b,lambda, u)=a^primeSigma_{12}b-frac{lambda}{2}(a^primeSigma_{11}a-1)-frac{ u}{2}(b^primeSigma_{22}b-1)label{equ:lagrange}end{equation}
对Lagrange函数 ef{equ:lagrange}求导得:
egin{equation}frac{partial L}{partial a}=Sigma_{12}b-lambdaSigma_{11}a=0label{equ:partiala}end{equation}
egin{equation}frac{partial L}{partial b}=Sigma_{21}a- uSigma_{22}b=0label{equ:partialb}end{equation}
将式子 ef{equ:partiala}左乘$a^prime$,式子 ef{equ:partialb}左乘$b^prime$得:
egin{align*}a^primeSigma_{12}b-lambda a^primeSigma_{11}a &=0\b^primeSigma_{21}a- u b^primeSigma_{22}b &=0end{align*}
又因为$(a^primeSigma_{12}b)^prime=b^primeSigma_{21}aLongrightarrow lambda a^primeSigma_{11}a= u b^primeSigma_{22}b$。由限制条件知:$lambda= u= ho=a^primeSigma_{12}b$,即$lambda$的值就是线性组合$U$和$V$的相关系数。我们重新将式子 ef{equ:partiala}和式子 ef{equ:partialb}写成:
egin{align}-lambdaSigma_{11}a+Sigma_{12}b&=0label{equ:partialaNew}\Sigma_{21}a-lambdaSigma_{22}b&=0label{equ:partialbNew}end{align}
然后将式子 ef{equ:partialbNew}左乘$Sigma_{12}Sigma_{22}^{-1}$得:
egin{equation}Sigma_{12}Sigma_{22}^{-1}Sigma_{21}a=lambdaSigma_{12}blabel{equ:partialbLeft}end{equation}
结合式子 ef{equ:partialaNew}得:
egin{equation}Sigma_{12}Sigma_{22}^{-1}Sigma_{21}a=lambda^2Sigma_{11}aLongrightarrow (Sigma_{12}Sigma_{22}^{-1}Sigma_{21}-lambda^2Sigma_{11})a=0label{equ:3to4}end{equation}
同理,将式子 ef{equ:partialaNew}左乘$Sigma_{21}Sigma_{11}^{-1}$,并将式子 ef{equ:partialbNew}代入式子 ef{equ:partialaNew}得:
egin{equation}(Sigma_{21}Sigma_{11}^{-1}Sigma_{12}-lambda^2Sigma_{22})b=0label{equ:4to3}end{equation}
将$Sigma_{11}^{-1}$左乘式子 ef{equ:3to4},$Sigma_{22}^{-1}$左乘式子 ef{equ:4to3}得:
egin{equation}left{egin{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=0end{array} ight. riangleqleft{egin{array}&mathbf{A}a=lambda^2a\mathbf{B}b=lambda^2bend{array} ight.end{equation}
说明:$lambda^2$既是矩阵$mathbf{A}$也是矩阵$mathbf{B}$的特征值,$a$与$b$分别是对应的特征向量。所以我们的问题转化成求矩阵$mathbf{A},mathbf{B}$的最大特征值对应的特征向量,而特征值的平方根$sqrt{lambda}$为相关系数,从而求出第一对典型相关变量。
此时,我们可以得到如下的猜想:是否矩阵$mathbf{A},mathbf{B}$的所有非零特征值的平方跟都会是其对应的典型相关系数?接下去,我们来证明如下猜想:
设$lambda^2_1geqlambda_2^2geqcdotsgeqlambda_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$对线性组合:
egin{equation}egin{array}&U_1=(a^{(1)})^primemathbf{X}&U_2=(a^{(2)})^primemathbf{X}&cdots&U_p=(a^{(p)})^primemathbf{X}\V_1=(b^{(1)})^primemathbf{Y}&V_2=(b^{(2)})^primemathbf{Y}&cdots&V_p=(b^{(p)})^primemathbf{Y}end{array}end{equation}
可以证明$(U_1,V_1),(U_2,V_2),cdots,(U_p,V_p)$就是其前$p$对典型变量,$lambda_1geqlambda_2geqcdotsgeqlambda_p$为其典型相关系数。
首先,在求出第一对典型变量的基础上求第二对典型变量。由上述分析我们可以知道该模型为:
egin{align}mathop{max}&quad (a^{(2)})^primeSigma_{12}b^{(2)} onumber\mathop{s.t.}&quad (a^{(2)})^primeSigma_{11} a^{(2)}=1 onumber\&quad(b^{(2)})^primeSigma_{22}b^{(2)}=1 onumber\&quad(a^{(1)})^primeSigma_{11}a^{(2)}=0label{con:1}\&quad(b^{(1)})^primeSigma_{22}b^{(2)}=0label{con:2}end{align}
其中限制 ef{con:1}和 ef{con:2}是由于第二对典型变量必须与第一对典型变量无关。其Lagrange方程以及相应的偏导为:
egin{align*}L(a^{(1)},b^{(1)},a^{(2)},b^{(2)}) &=(a^{(2)})^primeSigma_{12}b^{(2)}-frac{lambda}{2}((a^{(2)})^primeSigma_{11}a^{(2)}-1)-frac{ u}{2}((b^{(2)})^primeSigma_{22}b^{(2)}-1)\&quad + gamma (a^{(1)})^primeSigma_{11}a^{(2)}+eta(b^{(1)})^primeSigma_{22}b^{(2)}end{align*}
egin{align}frac{partial L}{partial a^{(2)}} &=Sigma_{12}b^{(2)}-lambdaSigma_{11}a^{(2)}+gammaSigma_{11}a^{(1)}=0label{partial:a2}\frac{partial L}{partial b^{(2)}}&=Sigma_{21}a^{(2)}- uSigma_{22}b^{(2)}+etaSigma_{22}b^{(1)}=0label{partial:b2}\frac{partial L}{partial a1}&=gammaSigma_{11}a^{(2)}=0label{partial:a1}\frac{partial L}{partial b^{(1)}}&=etaSigma_{22}b^{(2)}=0label{partial:b1}end{align}
将$(a^{(2)})^prime$左乘式子 ef{partial:a2},$(b^{(2)})^prime$左乘式子 ef{partial:b2}得:
egin{align}&(a^{(2)})^primeSigma_{12}b^{(2)}-lambda(a^{(2)})^primeSigma_{11}a^{(2)}+gamma (a^{(2)})^primeSigma_{11}a^{(1)}=0label{partial:a2new}\&(b^{(2)})^primeSigma_{21}a^{(2)}- u(b^{(2)})^primeSigma_{22}b^{(2)}+eta (b^{(2)})^primeSigma_{22}b^{(1)}=0label{partial:b2new}end{align}
将$(a^{(1)})^prime$左乘式子 ef{partial:a1},$(b^{(2)})^prime$左乘式子 ef{partial:b1}得:
egin{align}&(a^{(2)})^primeSigma_{11}a^{(1)}=0label{partial:a1new}\&(b^{(2)})^primeSigma_{22}b^{(2)}=0label{partial:b1new}end{align}
将式子 ef{partial:a1new}和式子 ef{partial:b1new}代入式子 ef{partial:a2new}和 ef{partial:b2new}得:
egin{align}&(a^{(2)})^primeSigma_{12}b^{(2)}-lambda(a^{(2)})^primeSigma_{11}a^{(2)}=0label{partial:a1a2}\&(b^{(2)})^primeSigma_{21}a^{(2)}- u(b^{(2)})^primeSigma_{22}b^{(2)}=0label{partial:b1b2}end{align}
其中式子 ef{partial:a1a2}和式子 ef{partial:b1b2}与式子 ef{equ:partiala}和式子 ef{equ:partialb}有相同的形式,只是$a,b$换成$a^{(2)},b^{(2)}$,故同样可以得到:
egin{equation}left{egin{array}&mathbf{A}a^{(2)}=lambda^2a^{(2)}\mathbf{B}b^{(2)}=lambda^2b^{(2)}end{array} ight.end{equation}
此时$a^{(2)} eq a,b^{(2)} eq b$,否则不满足限制 ef{con:1}和 ef{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$。