现象与现象直接的依存关系,从数据联系上看,可以分为两种不同的类型,即函数关系和相关关系。
函数关系是从数量上反映现象间的严格的依存关系,即与一个或几个变量取一定的值时,另一个变量有确定值与之相对应。相关关系是现象间不严格的依存关系,即各变量之间不存在确定性的关系。在相关关系中,当一个或几个相互联系的变量取一定数值时,与之相对应的另一变量值也发生相应变化,但其关系值不是固定的,往往按照某种规律在一定的范围内变化。
函数关系是指事物或现象之间存在着严格的依存关系,其主要特征是它的确定性,即对一个变量的每一个值,另一个变量都具有唯一确定的值与之相对应。相关关系反映出变量之间虽然相互影响,具有依存关系,但彼此之间是不能一一对应的。
如果两个变量间的相关程度不高,拟合回归方程便没有意义,因此相关分析往往在回归分析前进行。
1 基本介绍
典型相关分析(Canonical Correlation Analysis,以下简称CCA)是研究两组变量之间相关关系的一种多元统计方法。它能够揭示出两组变量之间的内在联系,是最常用的挖掘数据关联关系的算法之一。
在一元统计分析中,用相关系数来衡量两个随机变量之间的线性相关关系;用复相关系数研究一个随机变量和多个随机变量的线性相关关系。然而,这些统计方法在研究两组变量之间的相关关系时却无能为力。
在线性回归中,我们使用直线来拟合样本点,寻找
n
n
n 维特征向量
X
\pmb{X}
XXX 和输出结果(或者叫做label)
Y
Y
Y 之间的线性关系,其中
X
∈
R
n
,
Y
∈
R
\pmb{X}\in\mathbb{R}^n,Y\in\mathbb{R}
XXX∈Rn,Y∈R。然而当
Y
Y
Y 也是多维时,或者说
Y
Y
Y 也有多个特征时,我们希望分析出
X
\pmb{X}
XXX 和
Y
\pmb{Y}
YYY 的关系。
当然我们仍然可以使用回归的方法来分析,做法如下:
假设
X
∈
R
n
,
Y
∈
R
m
\pmb{X}\in\mathbb{R}^n,Y\in\mathbb{R}^m
XXX∈Rn,Y∈Rm,那么可以建立等式
Y
=
W
X
\pmb{Y}=\pmb{WX}
YYY=WXWXWX,如下
[ y 1 y 2 ⋮ y m ] = [ w 11 w 12 ⋯ w 1 n w 21 w 22 ⋯ w 2 n ⋮ ⋮ ⋮ ⋮ w m 1 w m 2 ⋯ w m n ] [ x 1 x 2 ⋮ x n ] (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)
其中 y i = w i T x y_i=\pmb{w}_i^T\pmb{x} yi=wwwiTxxx,形式和线性回归一样,需要训练 m m m 次得到 m m m 个 w i \pmb{w}_i wwwi。
这样做的一个缺点是, 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 各自特征间的线性组合,也就是考察 a T x \pmb{a}^T\pmb{x} aaaTxxx 和 b T y \pmb{b}^T\pmb{y} bbbTyyy 之间的关系。
在数理统计里面,我们都知道相关系数这个概念。假设有两个随机变量
X
X
X和
Y
Y
Y,则相关系数
ρ
\rho
ρ 的定义为:
ρ
(
X
,
Y
)
=
C
o
v
(
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)
其中
C
o
v
(
X
,
Y
)
Cov(X,Y)
Cov(X,Y) 是
X
X
X 和
Y
Y
Y 的协方差,而
μ
X
\mu_X
μX、
μ
Y
\mu_Y
μY分别是
X
X
X 和
Y
Y
Y 的均值,
D
(
X
)
D(X)
D(X)、
D
(
Y
)
D(Y)
D(Y) 分别是
X
X
X 和
Y
Y
Y 的方差。
相关系数
ρ
\rho
ρ 的取值为
[
−
1
,
1
]
[-1,1]
[−1,1],
ρ
\rho
ρ 的绝对值越接近于1,则
X
X
X 和
Y
Y
Y 的线性相关性越高;越接近于0,则
X
X
X 和
Y
Y
Y 的线性相关性越低。
典型相关分析的目的是识别并量化两组变量之间的联系,将两组变量相关关系的分析,转化为一组变量的线性组合与另一组变量线性组合之间的相关关系分析。
2 CCA 表示与求解
假设我们的数据集是
X
X
X 和
Y
Y
Y,
X
\pmb{X}
XXX 为
N
×
p
N\times p
N×p 的样本矩阵和
Y
\pmb{Y}
YYY 为
N
×
q
N\times q
N×q 的样本矩阵。其中
N
N
N 为样本个数,而
p
、
q
p、q
p、q 分别为
X
X
X 和
Y
Y
Y 的特征维度。此时,数据集为:
D
a
t
a
:
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
⋯
,
(
x
N
,
y
N
)
}
x
i
∈
R
p
,
y
i
∈
R
q
,
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
=
[
x
1
,
x
2
,
⋯
,
x
N
]
p
×
N
=
[
x
11
x
21
⋯
x
N
1
x
12
x
22
⋯
x
N
2
⋮
⋮
⋮
⋮
x
1
p
x
2
p
⋯
x
N
p
]
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 = [ y 1 , y 2 , ⋯ , y N ] q × N = [ y 11 y 21 ⋯ y N 1 y 12 y 22 ⋯ y N 2 ⋮ ⋮ ⋮ ⋮ y 1 q y 2 q ⋯ y N q ] 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 中有
N
N
N 个样本,每个样本
x
i
\pmb{x}_i
xxxi 为
p
p
p 维数据(含有
p
p
p 个feature
),数据
Y
\pmb{Y}
YYY 同样也有
N
N
N 个样本,每个样本
y
i
\pmb{y}_i
yyyi 为
q
q
q 维数据(含有
q
q
q 个feature
)。
对于
X
\pmb{X}
XXX 矩阵,我们将其投影到1维,或者说进行线性表示,对应的投影向量或者说线性系数向量为
a
i
∈
R
p
\pmb{a}_i \in \mathbb{R}^p
aaai∈Rp, 对于
Y
\pmb{Y}
YYY 矩阵,我们将其投影到1维,或者说进行线性表示,对应的投影向量或者说线性系数向量为
b
i
∈
R
q
\pmb{b}_i\in \mathbb{R}^q
bbbi∈Rq, 这样
X
,
Y
\pmb{X} ,\pmb{Y}
XXX,YYY 投影后得到的一维向量分别为
U
i
,
V
i
\pmb{U}_i,\pmb{V}_i
UUUi,VVVi(
i
=
1
,
⋯
,
m
i
n
{
p
,
q
}
i=1, \cdots, min\{p, q\}
i=1,⋯,min{p,q})。我们有
U
i
=
a
i
T
X
,
V
i
=
b
i
T
Y
(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
=
[
X
Y
]
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 自己的协方差矩阵;右上角是
C
o
v
(
X
,
Y
)
Cov(\pmb{X}, \pmb{Y})
Cov(XXX,YYY);左下角是
C
o
v
(
Y
,
X
)
Cov(\pmb{Y}, \pmb{X})
Cov(YYY,XXX),也是
Σ
12
\pmb{\Sigma}_{12}
ΣΣΣ12 的转置;右下角是
Y
\pmb{Y}
YYY 的协方差矩阵。
下面计算
U
i
\pmb{U}_i
UUUi 和
V
i
\pmb{V}_i
VVVi 的方差和协方差:
D
[
U
i
]
=
a
i
T
Σ
11
a
i
D
[
V
i
]
=
b
i
T
Σ
22
b
i
C
o
v
(
U
i
,
V
i
)
=
a
i
T
Σ
12
b
i
(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
[
U
i
]
=
D
[
a
i
T
X
]
=
1
N
∑
i
=
1
N
(
a
i
T
X
−
a
i
T
μ
X
)
2
=
a
i
T
1
N
∑
i
=
1
N
(
X
−
μ
X
)
2
a
i
=
a
i
T
Σ
11
a
i
(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)
下面就可以计算
ρ
(
U
i
,
V
i
)
\rho(\pmb{U}_i, \pmb{V}_i)
ρ(UUUi,VVVi)
ρ
(
U
i
,
V
i
)
=
a
i
T
Σ
12
b
i
a
i
T
Σ
11
a
i
b
i
T
Σ
22
b
i
(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)
(
U
i
,
V
i
)
(\pmb{U}_i, \pmb{V}_i)
(UUUi,VVVi) 为第
i
i
i 个canonical variate pair
,因为
p
≤
q
p\leq q
p≤q,所以就有
p
p
p 对canonical variate pair
。
我们希望找到使得在每一对canonical variate pair
中关联关系最大的线性组合。即:
a
r
g
m
a
x
⏟
a
,
b
ρ
(
U
i
,
V
i
)
=
a
r
g
m
a
x
⏟
a
,
b
a
i
T
Σ
12
b
i
a
i
T
Σ
11
a
i
b
i
T
Σ
22
b
i
(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,b
argmaxρ(UUUi,VVVi)=a,b
argmaxaaaiTΣΣΣ11aaaibbbiTΣΣΣ22bbbiaaaiTΣΣΣ12bbbi(2-9)
CCA 用于联合处理两组特征集,其背后的目标是寻找一对投影,每组各一个,使得在投影后得到的新特征最大程度地相关。把研究两组特征之间的问题化为研究两个所谓典型特征之间的相关问题。这里的典型特征不是从原特征组里挑出来的某个特征,而是原有特征的线性组合,因此需要求解的是这个线性组合的系数。
2.1 特征分解求解
在我们求导之前,需要对分母进行归一化,因为不做归一的话,
a
i
\pmb{a}_i
aaai 和
b
i
\pmb{b}_i
bbbi 扩大任何倍,上面的条件都成立,我们就无法确定
a
i
\pmb{a}_i
aaai 和
b
i
\pmb{b}_i
bbbi。因此我们令
a
i
T
Σ
11
a
i
=
b
i
T
Σ
22
b
i
=
1
\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i =\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i=1
aaaiTΣΣΣ11aaai=bbbiTΣΣΣ22bbbi=1(这里上下式子一同缩放不影响最后求
a
r
g
m
a
x
⏟
a
,
b
ρ
(
U
i
,
V
i
)
\underbrace{arg\;max}_{\boldsymbol{a},\boldsymbol{b}} \rho(\pmb{U}_i, \pmb{V}_i)
a,b
argmaxρ(UUUi,VVVi) 的结果)。此时优化问题的条件是:
M
a
x
i
m
i
z
e
:
a
i
T
Σ
12
b
i
S
u
b
j
e
c
t
t
o
:
a
i
T
Σ
11
a
i
=
1
,
b
i
T
Σ
22
b
i
=
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
)
=
a
i
T
Σ
12
b
i
−
λ
2
(
a
i
T
Σ
11
a
i
−
1
)
−
θ
2
(
b
i
T
Σ
22
b
i
−
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
=
Σ
12
b
i
−
λ
Σ
11
a
i
=
0
①
J
(
a
,
b
)
∂
b
=
Σ
21
a
i
−
θ
Σ
22
b
i
=
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)
等式①左乘
a
i
T
\pmb{a}_i^T
aaaiT,等式②左乘
b
i
T
\pmb{b}_i^T
bbbiT,再根据
a
i
T
Σ
11
a
i
=
b
i
T
Σ
22
b
i
=
1
\pmb{a}_i^T\pmb{\Sigma}_{11}\pmb{a}_i =\pmb{b}_i^T\pmb{\Sigma}_{22}\pmb{b}_i=1
aaaiTΣΣΣ11aaai=bbbiTΣΣΣ22bbbi=1,得到
λ
=
θ
=
a
i
T
Σ
12
b
i
(2-13)
\lambda=\theta=\pmb{a}_i^T\pmb{\Sigma}_{12}\pmb{b}_i \tag{2-13}
λ=θ=aaaiTΣΣΣ12bbbi(2-13)
其实也就是说我们的拉格朗日系数就是我们要优化的目标,所以只需要找到最大的
λ
\lambda
λ 即可。
我们上面的方程组进一步简化,并写成矩阵形式,得到
{
Σ
11
−
1
Σ
12
b
i
=
λ
a
i
③
Σ
22
−
1
Σ
21
a
i
=
λ
b
i
④
(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
−
1
0
0
Σ
22
−
1
]
[
0
Σ
12
Σ
21
0
]
[
a
i
b
i
]
=
λ
[
a
i
b
i
]
(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
=
[
Σ
11
0
0
Σ
22
]
,
A
=
[
0
Σ
12
Σ
21
0
]
,
w
=
[
a
i
b
i
]
(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
−
1
A
w
=
λ
w
(2-17)
\pmb{B}^{-1}\pmb{Aw}=\lambda\pmb{w} \tag{2-17}
BBB−1AwAwAw=λwww(2-17)
只要求得
B
−
1
A
\pmb{B}^{-1}\pmb{A}
BBB−1AAA 的最大特征值
λ
m
a
x
\lambda_{max}
λmax,那么
ρ
(
U
i
,
V
i
)
\rho(\pmb{U}_i, \pmb{V}_i)
ρ(UUUi,VVVi) 和
a
\pmb{a}
aaa 和
b
\pmb{b}
bbb 都可以求出。
如果直接计算
B
−
1
A
\pmb{B}^{-1}\pmb{A}
BBB−1AAA 的特征值,复杂度有点高。我们将式子④代入式子③,可以得到:
Σ
11
−
1
Σ
12
Σ
22
−
1
Σ
21
a
i
=
λ
2
a
i
(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 和特征向量 a i \pmb{a}_i aaai,然后根据式子④求得 b i \pmb{b}_i bbbi。
假设按照上述过程,得到了 λ \lambda λ 最大时的 a 1 \pmb{a}_1 aaa1 和 b 1 \pmb{b}_1 bbb1。那么 U 1 \pmb{U}_1 UUU1 和 V 1 \pmb{V}_1 VVV1 称为第一对典型变量(canonical variates), λ \lambda λ 即是 U 1 \pmb{U}_1 UUU1 和 V 1 \pmb{V}_1 VVV1 的相关系数。
最后,我们得到
U
1
\pmb{U}_1
UUU1 和
V
1
\pmb{V}_1
VVV1 的等式为:
U
1
=
a
1
T
X
,
V
1
=
b
1
T
Y
(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)
我们也可以接着去寻找第二组典型变量对,其最优化条件是
M
a
x
i
m
i
z
e
:
a
2
T
Σ
12
b
2
S
u
b
j
e
c
t
t
o
:
a
2
T
Σ
11
a
2
=
1
,
b
2
T
Σ
22
b
2
=
1
a
2
T
Σ
11
a
1
=
0
,
b
2
T
Σ
22
b
1
=
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)
第二组的约束条件等价于
ρ
(
U
2
,
U
1
)
=
0
,
ρ
(
V
2
,
V
1
)
=
0
\rho(\pmb{U}_2, \pmb{U}_1)=0,\rho(\pmb{V}_2, \pmb{V}_1)=0
ρ(UUU2,UUU1)=0,ρ(VVV2,VVV1)=0,主要是为了有效测度两组变量的相关信息,第二对典型变量应不再包含第一对典型变量已包含的信息。
计算步骤同第一组计算方法,只不过是
λ
2
2
\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 的第二大特征值。
类似地,依次可求出第 r r r 对典型变量: U r = a r T X \pmb{U}_r = \pmb{a}^T_r\pmb{X} UUUr=aaarTXXX 和 V r = b r T Y \pmb{V}_r = \pmb{b}^T_r\pmb{Y} VVVr=bbbrTYYY,其系数向量 a r \pmb{a}_r aaar 和 b r \pmb{b}_r bbbr 分别为矩阵 Σ 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 的第 r r r 特征根 λ r 2 \lambda_r^2 λr2 对应的特征向量。 λ r \lambda_r λr 即为第 r r 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 特征根及相应特征向量的求解。如果矩阵
Σ
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 的秩为
r
r
r ,则共有
r
r
r 对典型变量,第
k
k
k 对(
1
≤
k
≤
r
1\leq k \leq r
1≤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 第
k
k
k 特征根
λ
k
2
\lambda_k^2
λk2 相应的特征向量,典型相关系数为
λ
k
\lambda_k
λk。
典型变量具有如下性质:
- D ( U k ) = 1 , D ( V k ) = 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)
- C o v ( U i , U j ) = 0 , C o v ( V i , V j ) = 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)
- C o v ( U i , U j ) = { λ 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
/
2
u
,
b
=
Σ
22
−
1
/
2
v
\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
a
T
Σ
11
a
=
1
⇒
u
T
Σ
11
−
1
/
2
Σ
11
Σ
11
−
1
/
2
u
=
1
⇒
u
T
u
=
1
b
T
Σ
22
b
=
1
⇒
v
T
Σ
22
−
1
/
2
Σ
22
Σ
22
−
1
/
2
v
=
1
⇒
v
T
v
=
1
a
T
Σ
12
b
=
u
T
Σ
11
−
1
/
2
Σ
12
Σ
22
−
1
/
2
v
(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)
我们的优化目标变成下式:
a
r
g
m
a
x
⏟
u
,
v
u
T
Σ
11
−
1
/
2
Σ
12
Σ
22
−
1
/
2
v
S
u
b
j
e
c
t
t
o
:
u
T
u
=
1
,
v
T
v
=
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,v
argmaxuuuTΣΣΣ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
Σ
V
T
\pmb{M}=\pmb{U}\pmb{\Sigma}\pmb{V}^T
MMM=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 所有的列都为标准正交基,则
u
T
U
\pmb{u}^T\pmb{U}
uuuTUUU 和
V
T
v
\pmb{V}^T\pmb{v}
VVVTvvv 得到一个只有一个标量值为1,其余标量值为0的向量。此时我们有
u
T
Σ
11
−
1
/
2
Σ
12
Σ
22
−
1
/
2
v
=
u
T
U
Σ
V
T
v
=
σ
u
v
(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)
也就是说我们最大化 u T Σ 11 − 1 / 2 Σ 12 Σ 22 − 1 / 2 v \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 做了奇异值分解后,最大的奇异值就是我们优化目标的最大值,或者说我们的 U 1 \pmb{U}_1 UUU1 和 V 1 \pmb{V}_1 VVV1 之间的最大相关系数。利用对应的左右奇异向量 u \pmb{u} uuu、 v \pmb{v} vvv,我们也可以求出我们原始的 X \pmb{X} XXX 和 Y \pmb{Y} YYY 的线性系数 a = Σ 11 − 1 / 2 u , b = Σ 22 − 1 / 2 v \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()
![](https://img-blog.csdnimg.cn/f2c920d35c3841cdb686ea7391917716.png?x-oss-process=image/watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA6ZW_6Lev5ryr5ryrMjAyMQ==,size_20,color_FFFFFF,t_70,g_se,x_16#pic_center)
对于特征的意义,可以对照下图理解。
# 对数据标准化处理
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()
![](https://img-blog.csdnimg.cn/4055b220a7014e9eab5cde98e8c3956d.png?x-oss-process=image/watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA6ZW_6Lev5ryr5ryrMjAyMQ==,size_18,color_FFFFFF,t_70,g_se,x_16#pic_center)
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)
![](https://img-blog.csdnimg.cn/9cadf9deb42a4e719682500e43054d14.png?x-oss-process=image/watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA6ZW_6Lev5ryr5ryrMjAyMQ==,size_16,color_FFFFFF,t_70,g_se,x_16#pic_center)
从箱线图中可以清楚地看出,第一对典型协变量确实与种类变量高度相关。
最后,再将两个典型协变量一起绘制,以种类变量着色,将更清楚地揭示种类与两个典型协变量之间的关系。
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])
![](https://img-blog.csdnimg.cn/b51964ce01414594b15fcf193deb9ebe.png?x-oss-process=image/watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA6ZW_6Lev5ryr5ryrMjAyMQ==,size_16,color_FFFFFF,t_70,g_se,x_16#pic_center)
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
,
Z
X,Y,Z
X,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