典型相关分析 CCA

现象与现象直接的依存关系，从数据联系上看，可以分为两种不同的类型，即函数关系和相关关系。
函数关系是从数量上反映现象间的严格的依存关系，即与一个或几个变量取一定的值时，另一个变量有确定值与之相对应。相关关系是现象间不严格的依存关系，即各变量之间不存在确定性的关系。在相关关系中，当一个或几个相互联系的变量取一定数值时，与之相对应的另一变量值也发生相应变化，但其关系值不是固定的，往往按照某种规律在一定的范围内变化。

函数关系是指事物或现象之间存在着严格的依存关系，其主要特征是它的确定性，即对一个变量的每一个值，另一个变量都具有唯一确定的值与之相对应。相关关系反映出变量之间虽然相互影响，具有依存关系，但彼此之间是不能一一对应的。

如果两个变量间的相关程度不高，拟合回归方程便没有意义，因此相关分析往往在回归分析前进行。

1 基本介绍

典型相关分析（Canonical Correlation Analysis，以下简称CCA）是研究两组变量之间相关关系的一种多元统计方法。它能够揭示出两组变量之间的内在联系，是最常用的挖掘数据关联关系的算法之一。

在一元统计分析中，用相关系数来衡量两个随机变量之间的线性相关关系；用复相关系数研究一个随机变量和多个随机变量的线性相关关系。然而，这些统计方法在研究两组变量之间的相关关系时却无能为力。

在线性回归中，我们使用直线来拟合样本点，寻找 nnn 维特征向量 X\pmb{X}XXX 和输出结果（或者叫做label）YYY 之间的线性关系，其中 X∈Rn,Y∈R\pmb{X}\in\mathbb{R}^n,Y\in\mathbb{R}XXX∈Rn,Y∈R。然而当 YYY 也是多维时，或者说 YYY 也有多个特征时，我们希望分析出 X\pmb{X}XXX 和 Y\pmb{Y}YYY 的关系。
当然我们仍然可以使用回归的方法来分析，做法如下：
假设 X∈Rn,Y∈Rm\pmb{X}\in\mathbb{R}^n,Y\in\mathbb{R}^mXXX∈Rn,Y∈Rm，那么可以建立等式 Y=WX\pmb{Y}=\pmb{WX}YYY=WXWXWX，如下

[y1y2⋮ym]=[w11w12⋯w1nw21w22⋯w2n⋮⋮⋮⋮wm1wm2⋯wmn][x1x2⋮xn](1-1)\begin{bmatrix} y_1\\ y_2 \\ \vdots \\ y_m \end{bmatrix} = \begin{bmatrix} w_{11} & w_{12} & \cdots &w_{1n}\\ w_{21} & w_{22}& \cdots & w_{2n} \\ \vdots & \vdots & \vdots &\vdots \\ w_{m1} & w_{m2} &\cdots & w_{mn} \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \\ \vdots \\ x_n \end{bmatrix} \tag{1-1}⎣⎢⎢⎢⎡y1y2⋮ym⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡w11w21⋮wm1w12w22⋮wm2⋯⋯⋮⋯w1nw2n⋮wmn⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤(1-1)

其中 yi=wiTxy_i=\pmb{w}_i^T\pmb{x}yi=wwwiTxxx，形式和线性回归一样，需要训练 mmm 次得到 mmm 个 wi\pmb{w}_iwwwi。

这样做的一个缺点是，Y\pmb{Y}YYY 中的每个特征都与 X\pmb{X}XXX 的所有特征关联，Y\pmb{Y}YYY 中的特征之间没有什么联系。

我们想换一种思路来看这个问题，如果将 X\pmb{X}XXX 和 Y\pmb{Y}YYY 都看成整体，考察这两个整体之间的关系。我们将整体表示成 Y\pmb{Y}YYY 和 Y\pmb{Y}YYY 各自特征间的线性组合，也就是考察 aTx\pmb{a}^T\pmb{x}aaaTxxx 和 bTy\pmb{b}^T\pmb{y}bbbTyyy 之间的关系。

在数理统计里面，我们都知道相关系数这个概念。假设有两个随机变量 XXX和YYY，则相关系数 ρ\rhoρ 的定义为:
ρ(X,Y)=Cov(X,Y)D(X)D(Y)=E[(X−μX)(Y−μY)]D(X)D(Y)(1-2)\rho(X, Y)=\dfrac{Cov(X,Y)}{\sqrt{D(X)}\sqrt{D(Y)}}=\dfrac{E[(X-\mu_X)(Y-\mu_Y)]}{\sqrt{D(X)}\sqrt{D(Y)}} \tag{1-2}ρ(X,Y)=D(X)D(Y)Cov(X,Y)=D(X)D(Y)E[(X−μX)(Y−μY)](1-2)
其中 Cov(X,Y)Cov(X,Y)Cov(X,Y) 是 XXX 和 YYY 的协方差，而 μX\mu_XμX、μY\mu_YμY分别是 XXX 和 YYY 的均值，D(X)D(X)D(X)、D(Y)D(Y)D(Y) 分别是 XXX 和 YYY 的方差。
相关系数 ρ\rhoρ 的取值为 [−1,1][-1,1][−1,1],　ρ\rhoρ 的绝对值越接近于1，则 XXX 和 YYY 的线性相关性越高；越接近于0，则 XXX 和 YYY 的线性相关性越低。

典型相关分析的目的是识别并量化两组变量之间的联系，将两组变量相关关系的分析，转化为一组变量的线性组合与另一组变量线性组合之间的相关关系分析。

2 CCA 表示与求解

假设我们的数据集是 XXX 和 YYY， X\pmb{X}XXX 为 N×pN\times pN×p 的样本矩阵和 Y\pmb{Y}YYY 为 N×qN\times qN×q 的样本矩阵。其中 NNN 为样本个数，而 p、qp、qp、q 分别为 XXX 和 YYY 的特征维度。此时，数据集为：
Data：{(x1,y1),(x2,y2),⋯,(xN,yN)}xi∈Rp,yi∈Rq,i=1,2,⋯,N(2-1)Data：\{(\pmb{x}_1, \pmb{y}_1), (\pmb{x}_2, \pmb{y}_2), \cdots, (\pmb{x}_N, \pmb{y}_N)\} \quad \pmb{x}_{i}\in \mathbb{R}^{p},\pmb{y}_{i}\in \mathbb{R}^q,i=1,2,\cdots ,N \tag{2-1}Data：{(xxx1,yyy1),(xxx2,yyy2),⋯,(xxxN,yyyN)}xxxi∈Rp,yyyi∈Rq,i=1,2,⋯,N(2-1)
分别用矩阵表示：
X=[x1,x2,⋯,xN]p×N=[x11x21⋯xN1x12x22⋯xN2⋮⋮⋮⋮x1px2p⋯xNp]p×N(2-2)\pmb{X} = [\pmb{x}_{1},\pmb{x}_{2},\cdots ,\pmb{x}_{N}]_{p \times N} = \begin{bmatrix} x_{11} & x_{21} & \cdots & x_{N1} \\ x_{12} & x_{22} & \cdots & x_{N2} \\ \vdots & \vdots & \vdots &\vdots \\ x_{1p} & x_{2p} & \cdots & x_{Np} \\ \end{bmatrix}_{p \times N} \tag{2-2}XXX=[xxx1,xxx2,⋯,xxxN]p×N=⎣⎢⎢⎢⎡x11x12⋮x1px21x22⋮x2p⋯⋯⋮⋯xN1xN2⋮xNp⎦⎥⎥⎥⎤p×N(2-2)

Y=[y1,y2,⋯,yN]q×N=[y11y21⋯yN1y12y22⋯yN2⋮⋮⋮⋮y1qy2q⋯yNq]q×N(2-3)\pmb{Y} = [\pmb{y}_{1},\pmb{y}_{2},\cdots ,\pmb{y}_{N}]_{q \times N}= \begin{bmatrix} y_{11} & y_{21} & \cdots & y_{N1} \\ y_{12} & y_{22} & \cdots & y_{N2} \\ \vdots & \vdots & \vdots &\vdots \\ y_{1q} & y_{2q} & \cdots & y_{Nq} \\ \end{bmatrix}_{q \times N} \tag{2-3}YYY=[yyy1,yyy2,⋯,yyyN]q×N=⎣⎢⎢⎢⎡y11y12⋮y1qy21y22⋮y2q⋯⋯⋮⋯yN1yN2⋮yNq⎦⎥⎥⎥⎤q×N(2-3)

上面数据矩阵的解释：数据 X\pmb{X}XXX 中有 NNN 个样本，每个样本 xi\pmb{x}_ixxxi 为 ppp 维数据（含有 ppp 个feature），数据 Y\pmb{Y}YYY 同样也有 NNN 个样本，每个样本 yi\pmb{y}_iyyyi 为 qqq 维数据（含有 qqq 个feature）。

对于 X\pmb{X}XXX 矩阵，我们将其投影到1维，或者说进行线性表示，对应的投影向量或者说线性系数向量为 ai∈Rp\pmb{a}_i \in \mathbb{R}^paaai∈Rp, 对于 Y\pmb{Y}YYY 矩阵，我们将其投影到1维，或者说进行线性表示，对应的投影向量或者说线性系数向量为 bi∈Rq\pmb{b}_i\in \mathbb{R}^qbbbi∈Rq, 这样 X,Y\pmb{X} ,\pmb{Y}XXX,YYY 投影后得到的一维向量分别为 Ui,Vi\pmb{U}_i,\pmb{V}_iUUUi,VVVi（i=1,⋯,min{p,q}i=1, \cdots, min\{p, q\}i=1,⋯,min{p,q}）。我们有
Ui=aiTX,Vi=biTY(2-4)\pmb{U}_i = \pmb{a}_i^T\pmb{X}, \pmb{V}_i=\pmb{b}_i^T\pmb{Y} \tag{2-4}UUUi=aaaiTXXX,VVVi=bbbiTYYY(2-4)

形式化表示如下：
Z=[XY]E[Z]=[μXμY]Σ=D[Z]=[Σ11Σ12Σ21Σ22](2-5)\pmb{Z} = \begin{bmatrix} \pmb{X}\\ \pmb{Y}\end{bmatrix} \quad E[\pmb{Z}] = \begin{bmatrix} \pmb{\mu}_X\\ \pmb{\mu}_Y\end{bmatrix} \quad \pmb{\Sigma}=D[\pmb{Z}]= \begin{bmatrix} \pmb{\Sigma}_{11} & \pmb{\Sigma}_{12} \\ \pmb{\Sigma}_{21} & \pmb{\Sigma}_{22} \end{bmatrix} \tag{2-5}ZZZ=[XXXYYY]E[ZZZ]=[μμμXμμμY]ΣΣΣ=D[ZZZ]=[ΣΣΣ11ΣΣΣ21ΣΣΣ12ΣΣΣ22](2-5)

Σ\pmb{\Sigma}ΣΣΣ 是 Z\pmb{Z}ZZZ 的协方差矩阵；左上角是 X\pmb{X}XXX 自己的协方差矩阵；右上角是 Cov(X,Y)Cov(\pmb{X}, \pmb{Y})Cov(XXX,YYY)；左下角是 Cov(Y,X)Cov(\pmb{Y}, \pmb{X})Cov(YYY,XXX)，也是 Σ12\pmb{\Sigma}_{12}ΣΣΣ12 的转置；右下角是 Y\pmb{Y}YYY 的协方差矩阵。
下面计算 Ui\pmb{U}_iUUUi 和 Vi\pmb{V}_iVVVi 的方差和协方差：
D[Ui]=aiTΣ11aiD[Vi]=biTΣ22biCov(Ui,Vi)=aiTΣ12bi(2-6)D[\pmb{U}_i] = \pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i \quad D[\pmb{V}_i] = \pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i \quad Cov\left(\pmb{U}_i, \pmb{V}_i \right) = \pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i \tag{2-6}D[UUUi]=aaaiTΣΣΣ11aaaiD[VVVi]=bbbiTΣΣΣ22bbbiCov(UUUi,VVVi)=aaaiTΣΣΣ12bbbi(2-6)

推导一下第一个：
D[Ui]=D[aiTX]=1N∑i=1N(aiTX−aiTμX)2=aiT1N∑i=1N(X−μX)2ai=aiTΣ11ai(2-7)D[\pmb{U}_i]= D[\pmb{a}_i^T\pmb{X}] = \frac{1}{N}\sum_{i=1}^{N}(\pmb{a}_i^T\pmb{X}-\pmb{a}_i^T\pmb{\mu}_X)^2=\pmb{a}_i^T\frac{1}{N}\sum_{i=1}^{N}(\pmb{X}-\pmb{\mu}_X)^2\pmb{a}_i = \pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i \tag{2-7}D[UUUi]=D[aaaiTXXX]=N1i=1∑N(aaaiTXXX−aaaiTμμμX)2=aaaiTN1i=1∑N(XXX−μμμX)2aaai=aaaiTΣΣΣ11aaai(2-7)

下面就可以计算 ρ(Ui,Vi)\rho(\pmb{U}_i, \pmb{V}_i)ρ(UUUi,VVVi)
ρ(Ui,Vi)=aiTΣ12biaiTΣ11aibiTΣ22bi(2-8)\rho(\pmb{U}_i, \pmb{V}_i)=\dfrac{\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i}{\sqrt{\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i }\sqrt{\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i}} \tag{2-8}ρ(UUUi,VVVi)=aaaiTΣΣΣ11aaaibbbiTΣΣΣ22bbbiaaaiTΣΣΣ12bbbi(2-8)
(Ui,Vi)(\pmb{U}_i, \pmb{V}_i)(UUUi,VVVi) 为第 iii 个canonical variate pair，因为 p≤qp\leq qp≤q，所以就有 ppp 对canonical variate pair。
我们希望找到使得在每一对canonical variate pair中关联关系最大的线性组合。即：
argmax⏟a,bρ(Ui,Vi)=argmax⏟a,baiTΣ12biaiTΣ11aibiTΣ22bi(2-9)\underbrace{arg\;max}_{\boldsymbol{a},\boldsymbol{b}} \rho(\pmb{U}_i, \pmb{V}_i)=\underbrace{arg\;max}_{\boldsymbol{a},\boldsymbol{b}}\dfrac{\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i}{\sqrt{\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i }\sqrt{\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i}} \tag{2-9}a,bargmaxρ(UUUi,VVVi)=a,bargmaxaaaiTΣΣΣ11aaaibbbiTΣΣΣ22bbbiaaaiTΣΣΣ12bbbi(2-9)

CCA 用于联合处理两组特征集，其背后的目标是寻找一对投影，每组各一个，使得在投影后得到的新特征最大程度地相关。把研究两组特征之间的问题化为研究两个所谓典型特征之间的相关问题。这里的典型特征不是从原特征组里挑出来的某个特征，而是原有特征的线性组合，因此需要求解的是这个线性组合的系数。

2.1 特征分解求解

在我们求导之前，需要对分母进行归一化，因为不做归一的话，ai\pmb{a}_iaaai 和 bi\pmb{b}_ibbbi 扩大任何倍，上面的条件都成立，我们就无法确定 ai\pmb{a}_iaaai 和 bi\pmb{b}_ibbbi。因此我们令 aiTΣ11ai=biTΣ22bi=1\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i =\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i=1aaaiTΣΣΣ11aaai=bbbiTΣΣΣ22bbbi=1（这里上下式子一同缩放不影响最后求 argmax⏟a,bρ(Ui,Vi)\underbrace{arg\;max}_{\boldsymbol{a},\boldsymbol{b}} \rho(\pmb{U}_i, \pmb{V}_i)a,bargmaxρ(UUUi,VVVi) 的结果）。此时优化问题的条件是：
Maximize：aiTΣ12biSubjectto：aiTΣ11ai=1，biTΣ22bi=1(2-10)Maximize：\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i\\ Subject\quad to：\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i =1，\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i=1 \tag{2-10}Maximize：aaaiTΣΣΣ12bbbiSubjectto：aaaiTΣΣΣ11aaai=1，bbbiTΣΣΣ22bbbi=1(2-10)

这样使用拉格朗日乘子法就可以得到：
J(a,b)=aiTΣ12bi−λ2(aiTΣ11ai−1)−θ2(biTΣ22bi−1)(2-11)J(\pmb{a}, \pmb{b})=\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i - \frac{\lambda}{2}(\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i-1) - \frac{\theta}{2}(\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i-1) \tag{2-11}J(aaa,bbb)=aaaiTΣΣΣ12bbbi−2λ(aaaiTΣΣΣ11aaai−1)−2θ(bbbiTΣΣΣ22bbbi−1)(2-11)
        求偏导，并令导数为0，可得：
{J(a,b)∂a=Σ12bi−λΣ11ai=0①J(a,b)∂b=Σ21ai−θΣ22bi=0②(2-12)\begin{cases} \dfrac{J(\pmb{a}, \pmb{b})}{\partial \boldsymbol{a}} = \pmb{\Sigma}_{12}\pmb{b}_i - \lambda\pmb{\Sigma}_{11}\pmb{a}_i = 0 & \text {①}\\ \quad \\ \dfrac{J(\pmb{a}, \pmb{b})}{\partial \boldsymbol{b}} = \pmb{\Sigma}_{21}\pmb{a}_i - \theta\pmb{\Sigma}_{22}\pmb{b}_i=0 & \text {②}\end{cases} \tag{2-12}⎩⎪⎪⎪⎨⎪⎪⎪⎧∂aJ(aaa,bbb)=ΣΣΣ12bbbi−λΣΣΣ11aaai=0∂bJ(aaa,bbb)=ΣΣΣ21aaai−θΣΣΣ22bbbi=0①②(2-12)
        等式①左乘 aiT\pmb{a}_i^TaaaiT，等式②左乘 biT\pmb{b}_i^TbbbiT，再根据 aiTΣ11ai=biTΣ22bi=1\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i =\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i=1aaaiTΣΣΣ11aaai=bbbiTΣΣΣ22bbbi=1，得到
λ=θ=aiTΣ12bi(2-13)\lambda=\theta=\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i \tag{2-13}λ=θ=aaaiTΣΣΣ12bbbi(2-13)
        其实也就是说我们的拉格朗日系数就是我们要优化的目标，所以只需要找到最大的 λ\lambdaλ 即可。
        我们上面的方程组进一步简化，并写成矩阵形式，得到
　{Σ11−1Σ12bi=λai③Σ22−1Σ21ai=λbi④(2-14)\begin{cases}\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}\pmb{b}_i = \lambda\pmb{a}_i & \text { ③}\\ \quad \\ \pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{a}_i = \lambda\pmb{b}_i & \text {④}\end{cases} \tag{2-14}⎩⎪⎨⎪⎧ΣΣΣ11−1ΣΣΣ12bbbi=λaaaiΣΣΣ22−1ΣΣΣ21aaai=λbbbi ③④(2-14)

表示成矩阵形式
[Σ11−100Σ22−1][0Σ12Σ210][aibi]=λ[aibi](2-15)\begin{bmatrix} \pmb{\Sigma}_{11}^{-1} & 0 \\0 & \pmb{\Sigma}_{22}^{-1} \end{bmatrix}\begin{bmatrix} 0 & \pmb{\Sigma}_{12} \\ \pmb{\Sigma}_{21} & 0\end{bmatrix}\begin{bmatrix} \pmb{a}_i \\ \pmb{b}_i \end{bmatrix} = \lambda \begin{bmatrix} \pmb{a}_i \\ \pmb{b}_i \end{bmatrix} \tag{2-15}[ΣΣΣ11−100ΣΣΣ22−1][0ΣΣΣ21ΣΣΣ120][aaaibbbi]=λ[aaaibbbi](2-15)
        令
B=[Σ1100Σ22],A=[0Σ12Σ210],w=[aibi](2-16)\pmb{B}=\begin{bmatrix} \pmb{\Sigma}_{11}& 0 \\0 & \pmb{\Sigma}_{22} \end{bmatrix}, \pmb{A}=\begin{bmatrix} 0 & \pmb{\Sigma}_{12} \\ \pmb{\Sigma}_{21} & 0\end{bmatrix},\pmb{w} = \begin{bmatrix} \pmb{a}_i \\ \pmb{b}_i \end{bmatrix} \tag{2-16}BBB=[ΣΣΣ1100ΣΣΣ22],AAA=[0ΣΣΣ21ΣΣΣ120],www=[aaaibbbi](2-16)
        上式就可以写做：
B−1Aw=λw(2-17)\pmb{B}^{-1}\pmb{Aw}=\lambda\pmb{w} \tag{2-17}BBB−1AwAwAw=λwww(2-17)
        只要求得 B−1A\pmb{B}^{-1}\pmb{A}BBB−1AAA 的最大特征值 λmax\lambda_{max}λmax，那么 ρ(Ui,Vi)\rho(\pmb{U}_i, \pmb{V}_i)ρ(UUUi,VVVi) 和 a\pmb{a}aaa 和 b\pmb{b}bbb 都可以求出。

如果直接计算 B−1A\pmb{B}^{-1}\pmb{A}BBB−1AAA 的特征值，复杂度有点高。我们将式子④代入式子③，可以得到：
Σ11−1Σ12Σ22−1Σ21ai=λ2ai(2-18)\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12} \pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{a}_i =\lambda^2\pmb{a}_i \tag{2-18}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21aaai=λ2aaai(2-18)

这样先对 Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12} \pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 求特征值 λ2\lambda^2λ2 和特征向量 ai\pmb{a}_iaaai，然后根据式子④求得 bi\pmb{b}_ibbbi。

假设按照上述过程，得到了 λ\lambdaλ 最大时的 a1\pmb{a}_1aaa1 和 b1\pmb{b}_1bbb1。那么 U1\pmb{U}_1UUU1 和 V1\pmb{V}_1VVV1 称为第一对典型变量（canonical variates），λ\lambdaλ 即是 U1\pmb{U}_1UUU1 和 V1\pmb{V}_1VVV1 的相关系数。

最后，我们得到 U1\pmb{U}_1UUU1 和 V1\pmb{V}_1VVV1 的等式为：
U1=a1TX,V1=b1TY(2-19)\pmb{U}_1 = \pmb{a}_1^T\pmb{X}, \pmb{V}_1=\pmb{b}_1^T\pmb{Y} \tag{2-19}UUU1=aaa1TXXX,VVV1=bbb1TYYY(2-19)
我们也可以接着去寻找第二组典型变量对，其最优化条件是
Maximize：a2TΣ12b2Subjectto：a2TΣ11a2=1，b2TΣ22b2=1a2TΣ11a1=0，b2TΣ22b1=0(2-20)Maximize：\pmb{a}_2^T\pmb{\Sigma}_{12}\pmb{b}_2\\ Subject to：\pmb{a}_2^T\pmb{\Sigma}_{11}\pmb{a}_2 =1，\pmb{b}_2^T\pmb{\Sigma}_{22}\pmb{b}_2=1\\ \pmb{a}_2^T\pmb{\Sigma}_{11}\pmb{a}_1 =0，\pmb{b}_2^T\pmb{\Sigma}_{22}\pmb{b}_1=0 \tag{2-20}Maximize：aaa2TΣΣΣ12bbb2Subjectto：aaa2TΣΣΣ11aaa2=1，bbb2TΣΣΣ22bbb2=1aaa2TΣΣΣ11aaa1=0，bbb2TΣΣΣ22bbb1=0(2-20)

第二组的约束条件等价于 ρ(U2,U1)=0，ρ(V2,V1)=0\rho(\pmb{U}_2, \pmb{U}_1)=0，\rho(\pmb{V}_2, \pmb{V}_1)=0ρ(UUU2,UUU1)=0，ρ(VVV2,VVV1)=0，主要是为了有效测度两组变量的相关信息，第二对典型变量应不再包含第一对典型变量已包含的信息。
计算步骤同第一组计算方法，只不过是 λ22\lambda_2^2λ22 取 Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12} \pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 的第二大特征值。

类似地，依次可求出第 rrr 对典型变量：Ur=arTX\pmb{U}_r = \pmb{a}^T_r\pmb{X}UUUr=aaarTXXX 和 Vr=brTY\pmb{V}_r = \pmb{b}^T_r\pmb{Y}VVVr=bbbrTYYY，其系数向量 ar\pmb{a}_raaar 和 br\pmb{b}_rbbbr 分别为矩阵 Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 和 Σ22−1Σ21Σ11−1Σ12\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}ΣΣΣ22−1ΣΣΣ21ΣΣΣ11−1ΣΣΣ12 的第 rrr 特征根 λr2\lambda_r^2λr2 对应的特征向量。λr\lambda_rλr 即为第 rrr 典型相关系数。

综上所述，典型变量和典型相关系数的计算可归结为矩阵 Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 和 Σ22−1Σ21Σ11−1Σ12\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}ΣΣΣ22−1ΣΣΣ21ΣΣΣ11−1ΣΣΣ12 特征根及相应特征向量的求解。如果矩阵Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 和 Σ22−1Σ21Σ11−1Σ12\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}ΣΣΣ22−1ΣΣΣ21ΣΣΣ11−1ΣΣΣ12 的秩为 rrr ，则共有 rrr 对典型变量，第 kkk 对(1≤k≤r1\leq k \leq r1≤k≤r)典型变量的系数向量分别是矩阵Σ11−1Σ12Σ22−1Σ21\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}ΣΣΣ11−1ΣΣΣ12ΣΣΣ22−1ΣΣΣ21 和 Σ22−1Σ21Σ11−1Σ12\pmb{\Sigma}_{22}^{-1}\pmb{\Sigma}_{21}\pmb{\Sigma}_{11}^{-1}\pmb{\Sigma}_{12}ΣΣΣ22−1ΣΣΣ21ΣΣΣ11−1ΣΣΣ12 第 kkk 特征根λk2\lambda_k^2λk2 相应的特征向量，典型相关系数为 λk\lambda_kλk。
典型变量具有如下性质：

D(Uk)=1,D(Vk)=1(k=1,2,⋯,r)D(\pmb{U}_k)=1, D(\pmb{V}_k)=1(k=1, 2, \cdots, r)D(UUUk)=1,D(VVVk)=1(k=1,2,⋯,r)
Cov(Ui,Uj)=0,Cov(Vi,Vj)=0(i≠j)Cov(\pmb{U}_i, \pmb{U}_j)=0, Cov(\pmb{V}_i, \pmb{V}_j)=0(i\not=j)Cov(UUUi,UUUj)=0,Cov(VVVi,VVVj)=0(i=j)
Cov(Ui,Uj)={λi≠0(i=j,i=1,2,⋯,r)0(i≠j)Cov(\pmb{U}_i, \pmb{U}_j)= \begin{cases} \lambda_i\not=0\quad (i=j,i=1,2,\cdots,r) \\ 0\quad (i\not=j)\end{cases}Cov(UUUi,UUUj)={λi=0(i=j,i=1,2,⋯,r)0(i=j)

有些教程称以上方法为特征分解法，这里也可以使用奇异值分解（SVD）进行求解，下面写一下推导过程。

2.2 SVD求解

首先，令 a=Σ11−1/2u,b=Σ22−1/2v\pmb{a}=\pmb{\Sigma}_{11}^{-1/2}\pmb{u}, \pmb{b}=\pmb{\Sigma}_{22}^{-1/2}\pmb{v}aaa=ΣΣΣ11−1/2uuu,bbb=ΣΣΣ22−1/2vvv
aTΣ11a=1⇒uTΣ11−1/2Σ11Σ11−1/2u=1⇒uTu=1bTΣ22b=1⇒vTΣ22−1/2Σ22Σ22−1/2v=1⇒vTv=1aTΣ12b=uTΣ11−1/2Σ12Σ22−1/2v(2-21)\pmb{a}^T\pmb{\Sigma}_{11}\pmb{a} =1 \Rightarrow \pmb{u}^T\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{11}\pmb{\Sigma}_{11}^{-1/2}\pmb{u} =1 \Rightarrow \pmb{u}^T\pmb{u}=1\\ \quad \\ \pmb{b}^T\pmb{\Sigma}_{22}\pmb{b} =1 \Rightarrow \pmb{v}^T\pmb{\Sigma}_{22}^{-1/2}\pmb{\Sigma}_{22}\pmb{\Sigma}_{22}^{-1/2}\pmb{v} =1 \Rightarrow \pmb{v}^T\pmb{v}=1\\ \quad \\ \pmb{a}^T\pmb{\Sigma}_{12}\pmb{b} = \pmb{u}^T\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1/2}\pmb{v} \tag{2-21}aaaTΣΣΣ11aaa=1⇒uuuTΣΣΣ11−1/2ΣΣΣ11ΣΣΣ11−1/2uuu=1⇒uuuTuuu=1bbbTΣΣΣ22bbb=1⇒vvvTΣΣΣ22−1/2ΣΣΣ22ΣΣΣ22−1/2vvv=1⇒vvvTvvv=1aaaTΣΣΣ12bbb=uuuTΣΣΣ11−1/2ΣΣΣ12ΣΣΣ22−1/2vvv(2-21)
我们的优化目标变成下式：
argmax⏟u,vuTΣ11−1/2Σ12Σ22−1/2vSubjectto：uTu=1，vTv=1(2-22)\underbrace{arg\;max}_{\boldsymbol{u},\boldsymbol{v}}\quad\pmb{u}^T\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1/2}\pmb{v} \\ \quad \\ Subject\quad to：\pmb{u}^T\pmb{u}=1，\pmb{v}^T\pmb{v}=1\tag{2-22}u,vargmaxuuuTΣΣΣ11−1/2ΣΣΣ12ΣΣΣ22−1/2vvvSubjectto：uuuTuuu=1，vvvTvvv=1(2-22)
观察上面的式子，如果将 u\pmb{u}uuu 和 v\pmb{v}vvv 看做矩阵 M=Σ11−1/2Σ12Σ22−1/2\pmb{M}=\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1/2}MMM=ΣΣΣ11−1/2ΣΣΣ12ΣΣΣ22−1/2 的某一个奇异值对应的左右奇异向量。那么利用奇异值分解，我们可以得到 M=UΣVT\pmb{M}=\pmb{U}\pmb{\Sigma}\pmb{V}^TMMM=UUUΣΣΣVVVT，其中 U\pmb{U}UUU、V\pmb{V}VVV 分别为M\pmb{M}MMM 的左奇异向量和右奇异向量组成的矩阵，而 Σ\pmb{\Sigma}ΣΣΣ 为 M\pmb{M}MMM 的奇异值组成的对角矩阵。由于 U\pmb{U}UUU、V\pmb{V}VVV 所有的列都为标准正交基，则 uTU\pmb{u}^T\pmb{U}uuuTUUU 和 VTv\pmb{V}^T\pmb{v}VVVTvvv 得到一个只有一个标量值为1，其余标量值为0的向量。此时我们有
uTΣ11−1/2Σ12Σ22−1/2v=uTUΣVTv=σuv(2-23)\pmb{u}^T\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1/2}\pmb{v} = \pmb{u}^T\pmb{U\Sigma V}^T\pmb{v}=\sigma_{\boldsymbol{uv}}\tag{2-23}uuuTΣΣΣ11−1/2ΣΣΣ12ΣΣΣ22−1/2vvv=uuuTUΣVUΣVUΣVTvvv=σuv(2-23)

也就是说我们最大化 uTΣ11−1/2Σ12Σ22−1/2v\pmb{u}^T\pmb{\Sigma}_{11}^{-1/2}\pmb{\Sigma}_{12}\pmb{\Sigma}_{22}^{-1/2}\pmb{v}uuuTΣΣΣ11−1/2ΣΣΣ12ΣΣΣ22−1/2vvv，其实对应的最大值就是某一组左右奇异向量所对应的奇异值的最大值。也就是将 M\pmb{M}MMM 做了奇异值分解后，最大的奇异值就是我们优化目标的最大值，或者说我们的 U1\pmb{U}_1UUU1 和 V1\pmb{V}_1VVV1 之间的最大相关系数。利用对应的左右奇异向量 u\pmb{u}uuu、v\pmb{v}vvv，我们也可以求出我们原始的 X\pmb{X}XXX 和 Y\pmb{Y}YYY 的线性系数 a=Σ11−1/2u,b=Σ22−1/2v\pmb{a}=\pmb{\Sigma}_{11}^{-1/2}\pmb{u}, \pmb{b}=\pmb{\Sigma}_{22}^{-1/2}\pmb{v}aaa=ΣΣΣ11−1/2uuu,bbb=ΣΣΣ22−1/2vvv。

可以看出，SVD的求解方式非常简洁方便，但是两者求得的结果其实是等价的，只要利用SVD和特征分解之间的关系就很容易发现两者最后的结果相同。

3 Python案例

使用seaborn里提供的一个企鹅种类数据集来演示 CCA。

# 导入库文件
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.cross_decomposition import CCA

# 读取数据集
data = "https://raw.githubusercontent.com/mwaskom/seaborn-data/master/penguins.csv"
df = pd.read_csv(data)
df = df.dropna()   # Delete missing data in dataframe data
df.head()

对于特征的意义，可以对照下图理解。

# 对数据标准化处理
X = df[['bill_length_mm','bill_depth_mm']]
Y = df[['flipper_length_mm','body_mass_g']]
X_mc = (X-X.mean())/(X.std())
Y_mc = (Y-Y.mean())/(Y.std())

典型相关分析涉及多个变量，不同的变量往往具有不同的量纲及不同的数量级别。在进行典型相关分析时，由于典型变量是原始变量的线性组合，具有不同量纲变量的线性组合显然失去了实际意义。其次，不同的数量级别会导致“以大吃小”，即数量级别小的变量的影响会被忽略，从而影响了分析结果的合理性。因此，为了消除量纲和数量级别的影响，必须对数据先做标准化变换处理，然后再做典型相关分析。显然，经标准化变换之后的协差阵就是相关系数矩阵，因而，也即通常应从相关矩阵出发进行典型相关分析。

cca = CCA(n_components=2)
cca.fit(X_mc, Y_mc)
X_c, Y_c = cca.transform(X_mc, Y_mc)
scores = np.corrcoef(cca.x_scores_, cca.y_scores_, rowvar=False)
score = np.diag(scores[:2, 2:])
# score: array([0.78763151, 0.08638695])

现在我们已经完成了典型相关分析，让我们更深入地了解我们作为结果得到的典型协变量对。

在这个实验中，我们知道我们拥有的两组度量值，因为这两个数据矩阵来自同一组企鹅。我们早些时候怀疑这些测量的差异是由于企鹅物种差异造成的。因此，这两个测量值背后的一个共同潜在变量是物种变量。而我们的 CCA 分析的主要目标是捕捉共同变量。我们还看到第一对规范变量高度相关。

让我们检验一下典型协变量是否实际上是种类变量。首先，让我们用企鹅数据和第一对典型协变量创建数据矩。

# 用企鹅数据和第一对典型协变量创建数据矩阵
cc_res = pd.DataFrame({"CCX_1":X_c[:, 0],"CCY_1":Y_c[:, 0],"CCX_2":X_c[:, 1],"CCY_2":Y_c[:, 1],"Species":df.species.tolist(),"Island":df.island.tolist(),"sex":df.sex.tolist()})
cc_res.head()

plt.figure(figsize=(10, 6))
sns.boxplot(x="Species",y="CCX_1", data=cc_res)
sns.stripplot(x="Species",y="CCX_1", data=cc_res)

从箱线图中可以清楚地看出，第一对典型协变量确实与种类变量高度相关。

最后，再将两个典型协变量一起绘制，以种类变量着色，将更清楚地揭示种类与两个典型协变量之间的关系。

plt.figure(figsize=(10,6))
sns.scatterplot(x="CCX_1",y="CCY_1", hue="Species", data=cc_res)
plt.title('First Pair of Canonical Covariate, corr = %.2f' %np.corrcoef(X_c[:, 0], Y_c[:, 0])[0, 1])

4 小结

4.1 问题分析

若存在特征间线性相关时，协方差矩阵会存在不可逆的情况，该怎么处理？
一般遇到这种情况可以对 Σ11\pmb{\Sigma}_{11}ΣΣΣ11、Σ22\pmb{\Sigma}_{22}ΣΣΣ22 进行正则化，即变化为 Σ11+γI\pmb{\Sigma}_{11}+\gamma\pmb{I}ΣΣΣ11+γIII、Σ22+γI\pmb{\Sigma}_{22}+\gamma\pmb{I}ΣΣΣ22+γIII，然后继续求逆，其中 γ\gammaγ 为正则化系数。
样本量特别大时，该怎么办？
建议将样本对拆成两半的两个样本分别做典型相关分析，再把结果进行比较。
如果有多个集合(X,Y,ZX,Y,ZX,Y,Z)，怎么衡量多个样本集的关系？
可以使用文献里提到的典型相关的推广方法使得两两集合的距离差之和最小。

4.2 拓展

核典型相关分析(KCCA)
当我们的数据无法线性表示时，CCA就无法使用，此时我们可以利用核函数的思想，将数据映射到高维后，再利用CCA的思想降维到1维，求对应的相关系数和线性关系。
广义典型相关分析(GCCA)
GCCA通过改进CCA的优化条件，在内协方差矩阵中增加类别信息，使分类性能得以提高。
判别典型相关分析(DCCA)
DCCA则是将类别信息引入到互协方差矩阵，其充分考虑了同类样本之间的相关与不同类样本之间的相关的影响。
__

参考

典型相关分析 CCA 及 Python 实验：https://jishuin.proginn.com/p/763bfbd6a9ac
CCA典型关联分析原理与Python案例：https://cloud.tencent.com/developer/article/1652998
CCA算法-原理-python代码分析：https://ynqiu.github.io/2018/04/13/understanding-cca/
典型关联分析（Canonical Correlation Analysis）：https://www.cnblogs.com/jerrylead/archive/2011/06/20/2085491.html
典型相关分析：https://www.cnblogs.com/duye/p/9384821.html
典型关联分析（CCA）原理总结：https://www.cnblogs.com/pinard/p/6288716.html
GCCA、KCCA、TCCA、SCCA、DCCA：https://github.com/jameschapman19/cca_zoo