線性判別分析（Linear Discriminant Analysis）簡(jiǎn)稱LDA，是一種監(jiān)督學(xué)習(xí)方法。LDA是在目前機(jī)器學(xué)習(xí)、數(shù)據(jù)挖掘領(lǐng)域經(jīng)典且熱門(mén)的一個(gè)算法，據(jù)我所知，百度的商務(wù)搜索部里面就用了不少這方面的算法。
LDA的原理是，將帶上標(biāo)簽的數(shù)據(jù)（點(diǎn)），通過(guò)投影的方法，投影到維度更低的空間中，使得投影后的點(diǎn)，會(huì)形成按類別區(qū)分，一簇一簇的情況，相同類別的點(diǎn)，將會(huì)在投影后的空間中更接近。要說(shuō)明白LDA，首先得弄明白線性分類器：因?yàn)長(zhǎng)DA是一種線性分類器。對(duì)于K-分類的一個(gè)分類問(wèn)題，會(huì)有K個(gè)線性函數(shù)：
當(dāng)滿足條件：對(duì)于所有的j，都有Yk > Yj,的時(shí)候，我們就說(shuō)x屬于類別k。對(duì)于每一個(gè)分類，都有一個(gè)公式去算一個(gè)分值，在所有的公式得到的分值中，找一個(gè)最大的，就是所屬的分類了。
上式實(shí)際上就是一種投影，是將一個(gè)高維的點(diǎn)投影到一條高維的直線上，LDA最求的目標(biāo)是，給出一個(gè)標(biāo)注了類別的數(shù)據(jù)集，投影到了一條直線之后，能夠使得點(diǎn)盡量的按類別區(qū)分開(kāi)，當(dāng)k=2即二分類問(wèn)題的時(shí)候，如下圖所示：
clip_image002
紅色的方形的點(diǎn)為0類的原始點(diǎn)、藍(lán)色的方形點(diǎn)為1類的原始點(diǎn)，經(jīng)過(guò)原點(diǎn)的那條線就是投影的直線，從圖上可以清楚的看到，紅色的點(diǎn)和藍(lán)色的點(diǎn)被原點(diǎn)明顯的分開(kāi)了，這個(gè)數(shù)據(jù)只是隨便畫(huà)的，如果在高維的情況下，看起來(lái)會(huì)更好一點(diǎn)。下面我來(lái)推導(dǎo)一下二分類LDA問(wèn)題的公式：
假設(shè)用來(lái)區(qū)分二分類的直線（投影函數(shù))為：
image
LDA分類的一個(gè)目標(biāo)是使得不同類別之間的距離越遠(yuǎn)越好，同一類別之中的距離越近越好，所以我們需要定義幾個(gè)關(guān)鍵的值。
類別i的原始中心點(diǎn)為：（Di表示屬于類別i的點(diǎn))
image
類別i投影后的中心點(diǎn)為：
image
衡量類別i投影后，類別點(diǎn)之間的分散程度（方差）為：
image
最終我們可以得到一個(gè)下面的公式，表示LDA投影到w后的損失函數(shù)：
image
我們分類的目標(biāo)是，使得類別內(nèi)的點(diǎn)距離越近越好（集中），類別間的點(diǎn)越遠(yuǎn)越好。分 母表示每一個(gè)類別內(nèi)的方差之和，方差越大表示一個(gè)類別內(nèi)的點(diǎn)越分散，分子為兩個(gè)類別各自的中心點(diǎn)的距離的平方，我們最大化J(w)就可以求出最優(yōu)的w了。 想要求出最優(yōu)的w，可以使用拉格朗日乘子法，但是現(xiàn)在我們得到的J(w)里面，w是不能被單獨(dú)提出來(lái)的，我們就得想辦法將w單獨(dú)提出來(lái)。
我們定義一個(gè)投影前的各類別分散程度的矩陣，這個(gè)矩陣看起來(lái)有一點(diǎn)麻煩，其實(shí)意思是，如果某一個(gè)分類的輸入點(diǎn)集Di里面的點(diǎn)距離這個(gè)分類的中心店mi越近，則Si里面元素的值就越小，如果分類的點(diǎn)都緊緊地圍繞著mi，則Si里面的元素值越更接近0.
image
帶入Si，將J(w)分母化為：
image
image
同樣的將J(w)分子化為：
image
這樣損失函數(shù)可以化成下面的形式：
image
這樣就可以用最喜歡的拉格朗日乘子法了，但是還有一個(gè)問(wèn)題，如果分子、分母是都可以取任意值的，那就會(huì)使得有無(wú)窮解，我們將分母限制為長(zhǎng)度為1（這是用拉 格朗日乘子法一個(gè)很重要的技巧，在下面將說(shuō)的PCA里面也會(huì)用到，如果忘記了，請(qǐng)復(fù)習(xí)一下高數(shù)），并作為拉格朗日乘子法的限制條件，帶入得到：
image
這樣的式子就是一個(gè)求特征值的問(wèn)題了。
對(duì)于N(N>2)分類的問(wèn)題，我就直接寫(xiě)出下面的結(jié)論了：
image
這同樣是一個(gè)求特征值的問(wèn)題，我們求出的第i大的特征向量，就是對(duì)應(yīng)的Wi了。
這里想多談?wù)勌卣髦担卣髦翟诩償?shù)學(xué)、量子力學(xué)、固體力學(xué)、計(jì)算機(jī)等等領(lǐng)域都有廣泛的應(yīng)用，特征值表示的是矩陣的性質(zhì)，當(dāng)我們?nèi)〉骄仃嚨那癗個(gè)最大的特征 值的時(shí)候，我們可以說(shuō)提取到的矩陣主要的成分（這個(gè)和之后的PCA相關(guān)，但是不是完全一樣的概念）。在機(jī)器學(xué)習(xí)領(lǐng)域，不少的地方都要用到特征值的計(jì)算，比 如說(shuō)圖像識(shí)別、pagerank、LDA、還有之后將會(huì)提到的PCA等等。
下圖是圖像識(shí)別中廣泛用到的特征臉（eigen face），提取出特征臉有兩個(gè)目的，首先是為了壓縮數(shù)據(jù)，對(duì)于一張圖片，只需要保存其最重要的部分就是了，然后是為了使得程序更容易處理，在提取主要特 征的時(shí)候，很多的噪聲都被過(guò)濾掉了。跟下面將談到的PCA的作用非常相關(guān)。
image
特征值的求法有很多，求一個(gè)D * D的矩陣的時(shí)間復(fù)雜度是O(D^3), 也有一些求Top M的方法，比如說(shuō)power method，它的時(shí)間復(fù)雜度是O(D^2 * M), 總體來(lái)說(shuō)，求特征值是一個(gè)很費(fèi)時(shí)間的操作，如果是單機(jī)環(huán)境下，是很局限的。
#-*- coding: UTF-8 -*-
from numpy import *
import numpy
def lda(c1,c2):
#c1 第一類樣本，每行是一個(gè)樣本
#c2 第二類樣本，每行是一個(gè)樣本
#計(jì)算各類樣本的均值和所有樣本均值
m1=mean(c1,axis=0)#第一類樣本均值
m2=mean(c2,axis=0)#第二類樣本均值
c=vstack((c1,c2))#所有樣本
m=mean(c,axis=0)#所有樣本的均值
#計(jì)算類內(nèi)離散度矩陣Sw
n1=c1.shape[0]#第一類樣本數(shù)
print(n1);
n2=c2.shape[0]#第二類樣本數(shù)
#求第一類樣本的散列矩陣s1
s1=0
for i in range(0,n1):
s1=s1+(c1[i,:]-m1).T*(c1[i,:]-m1)
#求第二類樣本的散列矩陣s2
s2=0
for i in range(0,n2):
s2=s2+(c2[i,:]-m2).T*(c2[i,:]-m2)
Sw=(n1*s1+n2*s2)/(n1+n2)
#計(jì)算類間離散度矩陣Sb
Sb=(n1*(m-m1).T*(m-m1)+n2*(m-m2).T*(m-m2))/(n1+n2)
#求最大特征值對(duì)應(yīng)的特征向量
eigvalue,eigvector=linalg.eig(mat(Sw).I*Sb)#特征值和特征向量
indexVec=numpy.argsort(-eigvalue)#對(duì)eigvalue從大到小排序，返回索引
nLargestIndex=indexVec[:1] #取出最大的特征值的索引
W=eigvector[:,nLargestIndex] #取出最大的特征值對(duì)應(yīng)的特征向量
return W
function [ W] = FisherLDA(w1,w2)
%W最大特征值對(duì)應(yīng)的特征向量
%w1 第一類樣本
%w2 第二類樣本
%第一步：計(jì)算樣本均值向量
m1=mean(w1);%第一類樣本均值
m2=mean(w2);%第二類樣本均值
m=mean([w1;w2]);%總樣本均值
%第二步：計(jì)算類內(nèi)離散度矩陣Sw
n1=size(w1,1);%第一類樣本數(shù)
n2=size(w2,1);%第二類樣本數(shù)
%求第一類樣本的散列矩陣s1
s1=0;
for i=1:n1
s1=s1+(w1(i,:)-m1)'*(w1(i,:)-m1);
end
%求第二類樣本的散列矩陣s2
s2=0;
for i=1:n2
s2=s2+(w2(i,:)-m2)'*(w2(i,:)-m2);
end
Sw=(n1*s1+n2*s2)/(n1+n2);
%第三步：計(jì)算類間離散度矩陣Sb
Sb=(n1*(m-m1)'*(m-m1)+n2*(m-m2)'*(m-m2))/(n1+n2);
%第四步：求最大特征值和特征向量
%[V,D]=eig(inv(Sw)*Sb);%特征向量V，特征值D
A = repmat(0.1,[1,size(Sw,1)]);
B = diag(A);
[V,D]=eig(inv(Sw + B)*Sb);
[a,b]=max(max(D));
W=V(:,b);%最大特征值對(duì)應(yīng)的特征向量
end
cls1_data=[2.95 6.63;2.53 7.79;3.57 5.65;3.16 5.47];
cls2_data=[2.58 4.46;2.16 6.22;3.27 3.52];
%樣本投影前
plot(cls1_data(:,1),cls1_data(:,2),'.r');
hold on;
plot(cls2_data(:,1),cls2_data(:,2),'*b');
hold on;
W=FisherLDA(cls1_data,cls2_data);
%樣本投影后
new1=cls1_data*W;
new2=cls2_data*W;
k=W(2)/W(1);
plot([0,6],[0,6*k],'-k');
axis([2 6 0 11]);
hold on;
%畫(huà)出樣本投影到子空間點(diǎn)
for i=1:4
temp=cls1_data(i,:);
newx=(temp(1)+k*temp(2))/(k*k+1);
newy=k*newx;
plot(newx,newy,'*r');
end;
for i=1:3
temp=cls2_data(i,:);
newx=(temp(1)+k*temp(2))/(k*k+1);
newy=k*newx;
plot(newx,newy,'ob');
end;
