问题引入
现有三种鸢尾花分别标为0,1,2类,存在4种不同的属性,即萼片长(setosa),萼片宽(versicolor),花瓣长与花瓣的宽度,现在需要对这三种鸢尾花的分布在二维平面中进行可视化分析。数据来自于sklearn中的dataset部分数据。
萼片长 | 萼片宽 | 花瓣长 | 花瓣的宽度 | 类别 |
---|---|---|---|---|
5.1 | 3.5 | 1.4 | 0.2 | 0.0 |
4.9 | 3.0 | 1.4 | 0.2 | 0.0 |
4.7 | 3.2 | 1.3 | 0.2 | 0.0 |
7.0 | 3.2 | 4.7 | 1.4 | 1.0 |
6.4 | 3.2 | 4.5 | 1.5 | 1.0 |
6.9 | 3.1 | 4.9 | 1.5 | 1.0 |
6.3 | 3.3 | 6.0 | 2.5 | 2.0 |
5.8 | 2.7 | 5.1 | 1.9 | 2.0 |
7.1 | 3.0 | 5.9 | 2.1 | 2.0 |
6.3 | 2.9 | 5.6 | 1.8 | 2.0 |
高维数据进行可视化分析,常用思路是对数据进行降维。
符号说明
n
n
n表示列属性的个数。
x
i
\textbf{x}_i
xi为nx1的列向量,在表格中代表一个不包括类别属性的行。
c
=
0
,
1
,
⋯
,
k
,
⋯
c={0,1,\cdots,k,\cdots}
c=0,1,⋯,k,⋯表示类别的集合,对应上述表格中的0、1。
m
k
m^k
mk表示类别为c的样本总个数。
X
k
\textbf{X}^k
Xk表示类别k组成的矩阵,大小为
m
k
×
n
m^k\times{n}
mk×n。
x
i
k
,
k
∈
c
\textbf{x}_i^k,k\in{c}
xik,k∈c为nx1列向量,在表格中表示类别为k的属性行。
m
=
∑
k
∈
c
m
k
m=\sum_{k\in{c}}{m^k}
m=∑k∈cmk表示所有样本数。
θ
\theta
θ表示位于直线
l
l
l上的nx1的单位列向量。
前情提要
线性判别分析在二分类的应用中,我们有类内距离为
θ
T
S
w
θ
=
θ
T
∑
k
∈
c
∑
i
=
1
m
k
(
x
i
k
−
u
k
)
(
x
i
k
−
u
k
)
T
θ
\theta^T\textbf{S}_w\theta=\theta^T\sum_{k\in{c}}\sum_{i=1}^{m^k}(\textbf{x}_i^k-u^k)(\textbf{x}_i^k-u^k)^T\theta
θTSwθ=θTk∈c∑i=1∑mk(xik−uk)(xik−uk)Tθ
在看很多资料的时候,会发现人们习惯于直接定义类内散度为
S
w
=
∑
k
∈
c
∑
i
=
1
m
k
(
x
i
k
−
u
k
)
(
x
i
k
−
u
k
)
T
\textbf{S}_w=\sum_{k\in{c}}\sum_{i=1}^{m^k}(\textbf{x}_i^k-u^k)(\textbf{x}_i^k-u^k)^T
Sw=k∈c∑i=1∑mk(xik−uk)(xik−uk)T
然后再定义类内散度
S
b
=
(
u
0
−
u
1
)
(
u
0
−
u
1
)
T
\textbf{S}_b=(\textbf{u}^0-\textbf{u}^1)(\textbf{u}^0-\textbf{u}^1)^T
Sb=(u0−u1)(u0−u1)T
最后再在散度两边乘以直线向量
θ
\theta
θ得到映射后的类内距离与类间距离,这种定义初次看来毫无实际意义,但是熟悉PCA以及LDA的二分类之后,就发现这种定义方法所代表数学意义很有用,具体描述如下。
如果存在
Θ
=
(
θ
1
,
θ
2
,
.
.
.
,
θ
i
,
.
.
.
θ
n
)
\Theta=(\theta_1,\theta_2,...,\theta_i,...\theta_n)
Θ=(θ1,θ2,...,θi,...θn)组成一个新的标准正交基,那么
Θ
T
S
w
Θ
\Theta^T\textbf{S}_w\Theta
ΘTSwΘ就是新的标准正交基下的类内散度矩阵。对于对角线上的元素代表位于直线向量
θ
i
\theta_i
θi上的类内距离即
θ
i
S
w
θ
i
\theta_i\textbf{S}_w\theta_i
θiSwθi。
既然这种定义方式存在其意义,那么我们在本文中就采用这种直接定义散度矩阵的方式进行介绍。
推导过程
首先定义类内散度矩阵,与LDA的二分类一样
S
w
=
∑
k
∈
c
∑
i
=
1
m
k
(
x
i
k
−
u
k
)
(
x
i
k
−
u
k
)
T
(1)
\textbf{S}_w=\sum_{k\in{c}}\sum_{i=1}^{m^k}(\textbf{x}_i^k-u^k)(\textbf{x}_i^k-u^k)^T\tag1
Sw=k∈c∑i=1∑mk(xik−uk)(xik−uk)T(1)
再定义类间散度矩阵,与LDA稍微有些不同
S
b
=
∑
k
1
∈
c
∑
k
2
∈
c
m
k
1
m
k
2
(
u
k
1
−
u
k
2
)
(
u
k
1
−
u
k
2
)
T
=
∑
k
∈
c
m
k
(
u
k
−
u
)
(
u
k
−
u
)
T
(2)
\begin{matrix} \begin{aligned} \textbf{S}_b&=\sum_{k_1\in{c}}\sum_{k_2\in{c}}m^{k1}m^{k2}(\textbf{u}^{k_1}-\textbf{u}^{k_2})(\textbf{u}^{k_1}-\textbf{u}^{k_2})^T\\ &=\sum_{k\in{c}}m^k(\textbf{u}^k-\textbf{u})(\textbf{u}^k-\textbf{u})^T \end{aligned} \end{matrix}\tag2
Sb=k1∈c∑k2∈c∑mk1mk2(uk1−uk2)(uk1−uk2)T=k∈c∑mk(uk−u)(uk−u)T(2)
其中
u
=
1
m
∑
i
=
1
m
x
i
\textbf{u}=\frac{1}{m}\sum_{i=1}^{m}\textbf{x}_i
u=m1i=1∑mxi
每个平均向量可以代表该类,将每个类之间的散度矩阵相加,可以表示总体的散度矩阵,这些倒还可以理解,肯定会有人质疑前面为什么加上系数,为什么不定义成
S
b
=
∑
k
1
∈
c
∑
k
2
∈
c
(
u
k
1
−
u
k
2
)
(
u
k
1
−
u
k
2
)
T
\textbf{S}_b=\sum_{k_1\in{c}}\sum_{k_2\in{c}}(\textbf{u}^{k_1}-\textbf{u}^{k_2})(\textbf{u}^{k_1}-\textbf{u}^{k_2})^T
Sb=∑k1∈c∑k2∈c(uk1−uk2)(uk1−uk2)T形式,简单地阐述一下这种定义方式的合理性,我们知道二分类中定义成
S
b
=
(
u
0
−
u
1
)
(
u
0
−
u
1
)
T
\textbf{S}_b=(\textbf{u}^0-\textbf{u}^1)(\textbf{u}^0-\textbf{u}^1)^T
Sb=(u0−u1)(u0−u1)T,前面没有系数,
S
b
\textbf{S}_b
Sb就算前面加系数也只不过是对矩阵等比例放大而已,对结果并不会产生任何影响,但是对于类别较多的数据就不一样了,比如说0类与1类的训练数据有100个,而2类的训练数据只有1个,那么再求解的过程中将权重更加偏向于那些数据集本身较多的类别会变得更加合理,所以
S
b
\textbf{S}_b
Sb如此定义。
那么对于某个直线向量
θ
\theta
θ而言,类内距离为
θ
T
S
w
θ
(3)
\theta^T\textbf{S}_w\theta\tag3
θTSwθ(3)
类间距离为
θ
T
S
b
θ
(4)
\theta^T\textbf{S}_b\theta\tag4
θTSbθ(4)
求解目标仍然可以写为
arg
max
θ
θ
T
S
b
θ
θ
T
S
w
θ
(5)
\arg\max_{\theta}\frac{\theta^T\textbf{S}_b\theta}{\theta^T\textbf{S}_w\theta}\tag5
argθmaxθTSwθθTSbθ(5)
这种形式仍然可以用于分类,只不过是多分类罢了。
至于降维的应用,实际上还是需要数学理论的支持,首先明确一点,当数据从高维降到低维之后,各个维度之间一定是正交的,而且数据的特征会尽可能地被保留,LDA既然可以分类肯定可以保留特征,那么下面的推导就是证明LDA求解出来的多个极值解
θ
i
\theta_i
θi之间是正交的直线罢了。
具体内容见LDA二分类,构造拉格朗日求极值,其中
J
(
θ
)
=
θ
T
S
w
θ
+
λ
(
θ
T
S
b
θ
−
t
)
\textbf{J}(\theta)=\theta^T\textbf{S}_w\theta+\lambda(\theta^T\textbf{S}_b\theta-t)
J(θ)=θTSwθ+λ(θTSbθ−t)
对
θ
\theta
θ求导=0之后可以得到
S
b
θ
=
λ
S
w
θ
\textbf{S}_b\theta=\lambda\textbf{S}_{w}\theta
Sbθ=λSwθ
⟺
\iff
⟺
S
w
−
1
S
b
θ
=
λ
θ
(6)
\textbf{S}_{w}^{-1}\textbf{S}_b\theta=\lambda\theta\tag{6}
Sw−1Sbθ=λθ(6)
此时,
λ
\lambda
λ为
S
w
−
1
S
b
\textbf{S}_{w}^{-1}\textbf{S}_b
Sw−1Sb的特征值,
θ
\theta
θ为对应的特征向量,根据特征向量的性质,对于不同的
λ
i
\lambda_i
λi之间
θ
i
\theta_i
θi肯定是正交的,那么对于相同的
λ
\lambda
λ可以通过斯密特正交化方法正交化单位化处理成正交的形式,现在LDA可以进行降维才得到了数学上的保证。
对于多个
θ
\theta
θ而言,选取哪一个
θ
\theta
θ作为降维后的轴与PCA一样需要看
λ
\lambda
λ取值,证明如下:
S
w
\textbf{S}_w
Sw为实对称矩阵,那么存在
S
w
θ
=
t
1
θ
\textbf{S}_w\theta=t_1\theta
Swθ=t1θ
带入(6)中得到
S
b
θ
=
λ
t
1
θ
(7)
\textbf{S}_b\theta=\lambda{t_1}\theta\tag7
Sbθ=λt1θ(7)
与PCA选取较大的
λ
\lambda
λ一样,此处不再累述,在公式(6)中有,
λ
t
1
\lambda{t}_1
λt1越大,
θ
T
S
b
θ
\theta^T\textbf{S}_b\theta
θTSbθ越大,类间会越分散。
根据拉格朗日的约定条件
θ
T
S
w
θ
=
t
\theta^T\textbf{S}_w\theta=t
θTSwθ=t可以知道
t
1
t_1
t1为一个定值,且
t
1
=
t
t_1=t
t1=t。
所以,
λ
\lambda
λ越大,
θ
T
S
b
θ
\theta^T\textbf{S}_b\theta
θTSbθ越大,类间会越分散。
最终得到结论是选取
λ
\lambda
λ较大的矩阵。
那么对于降维而言,公式(5)就变成
arg
imax
θ
T
S
b
θ
θ
T
S
w
θ
(8)
\arg \textbf{imax}\space\frac{\theta^T\textbf{S}_b\theta}{\theta^T\textbf{S}_w\theta}\tag8
argimax θTSwθθTSbθ(8)
其中imax代表求极值
求解方式与公式(6)一样,如下
S
w
−
1
S
b
θ
=
λ
θ
\textbf{S}_{w}^{-1}\textbf{S}_b\theta=\lambda\theta
Sw−1Sbθ=λθ
如果需要降维到d维,只需要求出前d个较大的
λ
\lambda
λ所对应的
θ
\theta
θ即可。
补充知识
1. 为什么求解公式定义成
arg
max
θ
θ
T
S
b
θ
θ
T
S
w
θ
\arg\max_{\theta}\frac{\theta^T\textbf{S}_b\theta}{\theta^T\textbf{S}_w\theta}
argθmaxθTSwθθTSbθ
却不直接定义成
arg
min
θ
θ
T
S
w
θ
θ
T
S
b
θ
\arg\min_{\theta}\frac{\theta^T\textbf{S}_w\theta}{\theta^T\textbf{S}_b\theta}
argθminθTSbθθTSwθ
实际上,求最小值的方式也对,利用拉格朗日求解之后就会发现,
λ
\lambda
λ越小,
θ
T
S
w
θ
\theta^T\textbf{S}_w\theta
θTSwθ越小,我们目标是求解类内距离最小,最终降维后选取的是前d个较小的
λ
\lambda
λ,与PCA结论恰好相反,可能为了方便记忆,故意写成第一个公式进行求解。
2. 公式(8)是本人根据自己理解写出的公式,实际上教材中的公式常常写作
arg
max
t
r
(
Θ
T
S
b
Θ
)
t
r
(
Θ
T
S
w
Θ
)
\arg\max\frac{tr(\Theta^T\textbf{S}_b\Theta)}{tr(\Theta^T\textbf{S}_w\Theta)}
argmaxtr(ΘTSwΘ)tr(ΘTSbΘ)
Θ
=
(
θ
1
,
θ
2
,
.
.
.
,
θ
i
,
.
.
.
θ
n
)
\Theta=(\theta_1,\theta_2,...,\theta_i,...\theta_n)
Θ=(θ1,θ2,...,θi,...θn)
这种写法比较简洁,但是对于刚刚入门的人而言可能有点难以理解,求迹就是求对角线和,等价于
arg
max
∑
i
=
1
n
(
θ
i
T
S
b
θ
i
)
∑
i
=
1
n
(
θ
i
T
S
w
θ
i
)
\arg\max\frac{\sum_{i=1}^n(\theta_i^T\textbf{S}_b\theta_i)}{\sum_{i=1}^n(\theta_i^T\textbf{S}_w\theta_i)}
argmax∑i=1n(θiTSwθi)∑i=1n(θiTSbθi)
写成拉格朗日的形式就会发现,对
θ
i
\theta_i
θi求导,那么
θ
j
,
i
≠
j
\theta_j,i\neq{j}
θj,i=j相当于是个常数,求导之后含有
θ
j
\theta_j
θj的项变成0,结果不变。
3. 构造拉格朗日的时候,本人将分母写作常数t,很多教材都写作1,求导结果都一样,等于1的写法,是数学中的一个求解定理,需要拥有相应的数学知识才能立即反应过来,实际上
θ
\theta
θ是个直线,也就是说,不管是
θ
\theta
θ还是
t
θ
\sqrt{t}\theta
tθ都是正确解,所以分母等于常数t还是等于1,都可以约掉,本人希望尽可能将数学知识限定在非数学系大学本科数学教学范围之内,不愿做过多的拓展,故意写作成常数t。
python代码
利用sklearn降维之后的图如下,代码来自于sklearn官方文档
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
iris = datasets.load_iris()
X = iris.data
y = iris.target
target_names = iris.target_names
lda = LDA(n_components=2)
X_r2 = lda.fit(X, y).transform(X)
plt.figure()
for c, i, target_name in zip("rgb", [0, 1, 2], target_names):
plt.scatter(X_r2[y == i, 0], X_r2[y == i, 1], c=c, label=target_name)
plt.legend()
plt.title('LDA of IRIS dataset')
plt.show()
下面是自己的python代码实现结果,虽然形式上有区别,但两者都可以区分三个不同种类
源码如下
from sklearn import datasets
import numpy as np
import matplotlib.pyplot as plt
#有num个种类
#降到k维
def desc(num,k):
iris = datasets.load_iris()
xa = iris.data
y = iris.target
target_names = iris.target_names
#所有元素的平均值
ua=np.array([np.mean(xa,axis=0)])
n=xa.shape[1]
#存储x值
x=[[]for i in range(num)]
u=[[]for i in range(num)]
sw=np.zeros([n,n])
sb=np.zeros([n,n])
for i in range(y.shape[0]):
x[y[i]].append(xa[i])
for i in range(num):
x[i]=np.array(x[i])
u[i]=np.array([np.mean(x[i],axis=0)])
#去中心化,计算Sw
for i in range(num):
x[i]=x[i]-u[i]
sw=sw+np.dot(x[i].T,x[i])
sb=sb+x[i].shape[0]*np.dot((u[i]-ua).T,(u[i]-ua))
#求Sw_1Sb的特征值与特征向量
eig_value,eig_vec=np.linalg.eig(np.dot(np.linalg.inv(sw),sb))
print(eig_vec)
#对特征值进行排序,返回标号
index=np.argsort(eig_value)
#获取取最后k个特征向量
eig_vec_k=eig_vec[:,index[:-k-1:-1]]
colors=["r","g","b"]
#绘图
for i in range(xa.shape[0]):
plt.scatter(np.dot(xa[i],eig_vec_k[:,0]),np.dot(xa[i],eig_vec_k[:,1]),color=colors[y[i]])
plt.title("my LDA")
plt.show()
desc(3,2)