3-1 二类LDA算法及MATLAB实现
1.二类LDA原理
\quad \quad
LDA(Linear Discriminant Analysis)是是一种监督学习的降维技术。LDA的思想可以用一句话概括,就是“投影后类内方差最小,类间方差最大”,即分类依据是(协)方差。
注 在概率论和统计学中,协方差用于衡量两个变量的总体误差,设
x
,
y
x,y
x,y是两个随机变量,则协方差定义为:
c
o
v
(
x
,
y
)
=
E
[
(
x
−
E
(
x
)
)
(
y
−
E
(
y
)
)
]
\quad \quad\quad\quad\quad cov(x,y)=E[(x-E(x))(y-E(y))]
cov(x,y)=E[(x−E(x))(y−E(y))]
\quad \quad
对于二类LDA而言,它的输出结果是二维的。我们希望给定的数据投影到一维的一条直线,让每一种类别数据的投影点尽可能的接近。
\quad \quad
假设数据集
D
=
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
m
,
y
m
)
D={(x_1,y_1),(x_2,y_2),…,(x_m,y_m)}
D=(x1,y1),(x2,y2),…,(xm,ym),其中任意样本
x
i
x_i
xi为
n
n
n维向量,
y
i
∈
{
0
,
1
}
y_i\in\{0,1\}
yi∈{0,1}。我们定义
N
j
,
j
=
0
,
1
N_j,j=0,1
Nj,j=0,1为第
j
j
j类样本的个数(两数值一般相同),
X
j
X_j
Xj为第
j
j
j类样本的集合,而
u
j
(
j
=
0
,
1
)
u_j(j=0,1)
uj(j=0,1)为第
j
j
j类样本的均值向量,定义
∑
j
(
j
=
0
,
1
)
\sum_j(j=0,1)
∑j(j=0,1)为第
j
j
j类样本的协方差矩阵(严格说是缺少分母部分的协方差矩阵).则:
\quad \quad
\quad \quad
u
j
=
1
N
j
∑
x
∈
X
j
x
u_j=\frac{1}{N_j}\underset{x\in X_j}{\sum}x
uj=Nj1x∈Xj∑x
\quad \quad
\quad \quad
∑
j
=
1
N
j
∑
x
∈
X
j
(
x
−
u
j
)
(
x
−
u
j
)
T
\sum_j=\frac{1}{N_j}\underset{x\in X_j}{\sum}(x-u_j)(x-u_j)^T
∑j=Nj1x∈Xj∑(x−uj)(x−uj)T
由于输出结果只有两类,我们将数据投影到一条直线上即可。
于是,分类的关键在于确定投影直线。不妨设投影直线是向量
w
w
w,则对任意一个样本
x
i
x_i
xi,它在直线
w
w
w的投影为
w
T
x
i
w^Tx_i
wTxi,对于两个类别的中心点
u
0
u_0
u0,
u
1
u_1
u1,在直线
w
w
w的投影为
w
T
u
0
w^Tu_0
wTu0和
w
T
u
1
w^Tu_1
wTu1。若要使得类间距离最大,即使得:
\quad \quad\quad \quad
∣
∣
w
T
u
0
−
w
T
u
1
∣
∣
2
2
||w^Tu_0-w^Tu_1||^{2}_{2}
∣∣wTu0−wTu1∣∣22
最大化。若要使得同一种类别数据的投影点尽可能的接近,也就是要同类样本投影点的协方差:
\quad \quad\quad \quad
w
T
∑
0
w
w^T\sum_0w
wT∑0w,
\quad
w
T
∑
1
w
w^T\sum_1w
wT∑1w
尽可能的小,即最小化
\quad \quad\quad \quad
w
T
∑
0
w
+
w
T
∑
0
w
w^T\sum_0w+w^T\sum_0w
wT∑0w+wT∑0w。
于是,优化目标变为:
\quad
a
r
g
m
a
x
J
(
w
)
=
∣
∣
w
T
u
0
−
w
T
u
1
∣
∣
2
2
w
T
∑
0
w
+
w
T
∑
0
w
=
w
T
(
u
0
−
u
1
)
(
u
0
−
u
1
)
T
w
w
T
∑
0
w
+
w
T
∑
0
w
argmaxJ(w)=\frac{||w^Tu_0-w^Tu_1||^{2}_{2}}{w^T\sum_0w+w^T\sum_0w}=\frac{w^T(u_0-u_1)(u_0-u_1)^Tw}{w^T\sum_0w+w^T\sum_0w}
argmaxJ(w)=wT∑0w+wT∑0w∣∣wTu0−wTu1∣∣22=wT∑0w+wT∑0wwT(u0−u1)(u0−u1)Tw
为了计算优化目标,引入瑞利商的定义:
瑞利商是形如
\quad \quad\quad \quad
R
(
A
,
x
)
=
x
T
A
x
x
T
x
R(A,x)=\frac{x^TAx}{x^Tx}
R(A,x)=xTxxTAx
的函数,其中
x
x
x是n维向量,
A
A
A是
n
×
n
n\times n
n×n的Hermitan矩阵,则瑞利商满足即它的最大值等于矩阵
A
A
A最大的特征值,而最小值等于矩阵
A
A
A的最小的特征值。
证:由正交变换的知识, x T A x x^TAx xTAx可相似对角化为 λ 1 y 1 2 + . . . + λ n y n 2 \lambda_1y_1^2+...+\lambda_ny_n^2 λ1y12+...+λnyn2,则 λ m i n ( y 1 2 + . . . + y n 2 ) ≤ λ 1 y 1 2 + . . . + λ n y n 2 ≤ λ m a x ( y 1 2 + . . . + y n 2 ) \lambda_{min}(y_1^2+...+y_n^2)\leq\lambda_1y_1^2+...+\lambda_ny_n^2\leq\lambda_{max}(y_1^2+...+y_n^2) λmin(y12+...+yn2)≤λ1y12+...+λnyn2≤λmax(y12+...+yn2)
回到本文,令类内散度矩阵
S
w
S_w
Sw为
S
w
=
∑
0
+
∑
1
=
∑
x
∈
X
0
(
x
−
x
0
)
(
x
−
x
0
)
T
+
∑
x
∈
X
1
(
x
−
x
1
)
(
x
−
x
1
)
T
S_w=\sum_0+\sum_1=\underset{x\in X_0}{\sum}(x-x_0)(x-x_0)^T+\underset{x\in X_1}{\sum}(x-x_1)(x-x_1)^T
Sw=∑0+∑1=x∈X0∑(x−x0)(x−x0)T+x∈X1∑(x−x1)(x−x1)T
同时定义类间散度矩阵
S
b
S_b
Sb为
\quad \quad\quad \quad
S
b
=
(
u
0
−
u
1
)
(
u
0
−
u
1
)
T
S_b=(u_0-u_1)(u_0-u_1)^T
Sb=(u0−u1)(u0−u1)T
则优化目标变为
\quad \quad\quad \quad
a
r
g
m
a
x
J
(
x
)
=
w
T
S
b
w
w
T
S
w
w
argmaxJ(x)=\frac{w^TS_bw}{w^TS_{w}w}
argmaxJ(x)=wTSwwwTSbw
显然,
S
w
S_w
Sw正定,利用瑞利商的性质,
J
(
x
)
J(x)
J(x)的最大值为
S
w
−
1
2
S
b
S
w
−
1
2
S_w^{-\frac{1}{2}}S_bS_w^{-\frac{1}{2}}
Sw−21SbSw−21的最大特征值。对应的
w
w
w是
S
w
−
1
2
S
b
S
w
−
1
2
S_w^{-\frac{1}{2}}S_bS_w^{-\frac{1}{2}}
Sw−21SbSw−21的最大特征值所对应的特征向量。
此时
S
w
−
1
S
b
S_w^{-1}S_b
Sw−1Sb的特征值与
S
w
−
1
2
S
b
S
w
−
1
2
S_w^{-\frac{1}{2}}S_bS_w^{-\frac{1}{2}}
Sw−21SbSw−21的特征值相等,
S
w
−
1
S
b
S_w^{-1}S_b
Sw−1Sb的最大特征值所对应的特征向量
w
′
=
S
w
−
1
2
w
w'=S_w^{-\frac{1}{2}}w
w′=Sw−21w。
证明:设
A
=
S
w
−
1
2
S
b
S
w
−
1
2
A=S_w^{-\frac{1}{2}}S_bS_w^{-\frac{1}{2}}
A=Sw−21SbSw−21,
B
=
S
w
−
1
S
b
B=S_w^{-1}S_b
B=Sw−1Sb,设
λ
,
x
\lambda,x
λ,x分别为
A
A
A的特征值和特征向量,则
A
x
=
S
w
−
1
2
S
b
S
w
−
1
2
x
=
λ
x
\quad \quad\quad\quad Ax = S_w^{-\frac{1}{2}}S_bS_w^{-\frac{1}{2}}x=\lambda x
Ax=Sw−21SbSw−21x=λx
等式两端同时左乘
S
w
−
1
2
S_w^{-\frac{1}{2}}
Sw−21有:
S
w
−
1
2
A
x
=
S
w
−
1
S
b
S
w
−
1
2
x
=
S
w
−
1
2
λ
x
=
λ
S
w
−
1
2
x
\quad \quad\quad S_w^{-\frac{1}{2}}Ax = S_w^{-1}S_bS_w^{-\frac{1}{2}}x=S_w^{-\frac{1}{2}}\lambda x=\lambda S_w^{-\frac{1}{2}}x
Sw−21Ax=Sw−1SbSw−21x=Sw−21λx=λSw−21x
⟺
\iff
⟺
B
S
w
−
1
2
x
=
λ
S
w
−
1
2
x
\quad \quad\quad BS_w^{-\frac{1}{2}}x =\lambda S_w^{-\frac{1}{2}}x
BSw−21x=λSw−21x
从而
λ
,
S
w
−
1
2
x
\lambda,S_w^{-\frac{1}{2}}x
λ,Sw−21x分别是
B
B
B的特征值和特征向量。
2.算例与程序
其他文章中讲的LDA算法实例也都不错,但由于不太好找例子,我只好用MATLAB及Excel造一个例子。
给定一个随机向量
A
A
A,对随机生成的其他向量,若
A
x
>
c
Ax>c
Ax>c(c为任意常数),则类别分为1,反之分为0.显然分类依据是线形的。
令A=0,2,3,-4
B=
加上分类后的数据
编写程序求投影向量 w w w
clc,clear
A=xlsread('test1.xlsx');
n=0;m=0;
[s,t]=size(A(:,1:end-1));
sum1=zeros(1,t);sum2=sum1;
for i=1:s
if A(i,end)==1
n=n+1;
sum1=sum1+A(i,1:end-1);
else
m=m+1;
sum2=sum2+A(i,1:end-1);
end
end
u1=sum1./n;
u2=sum2./m;
u1=u1';u2=u2';
a1=zeros(t,t);a2=a1;
for i=1:s
if A(i,end)==1
a1=a1+(A(i,1:end-1)'-u1)*(A(i,1:end-1)'-u1)';
else
a2=a2+(A(i,1:end-1)'-u2)*(A(i,1:end-1)'-u2)';
end
end
sw=a1+a2;
sb=(u1-u2)*(u1-u2)';
[x,y]=eig(sw\sb);
eigenvalue=diag(y);%求对角线向量
lamda=max(eigenvalue);%求最大特征值
for i=1:length(y)%求最大特征值对应的序数
if lamda==eigenvalue(i)
break;
end
end
[m,n]=size(x);%得到行数和列数
y_lamda=x(:,i);%求矩阵最大特征值对应的特征向量
w=sw^(-1/2)*y_lamda;
验证:
随机生成一个向量
x
1
x1
x1
x1=rand(10,1);
x=
用A作用x的判别结果(结果大于1判为1,否则为0)为(0)
用w作用x的判别结果
用w作用后为原矩阵(B)被将为1维
我们发现正值在原判别中是0,负值在原判别中是1,随机生成的x被w作用的值为正,判别结果应该为0,与被A作用结果相同。