机器学习算法之LDA(线性判别分析)人脸识别FisherFaces,最大化类间散度,最小化类内散度,LDA线性判别分析原理剖析,LDA人脸识别matlab实现,LDA人脸识别python实现

LDA算法概念

       线性判别分析(Linear discriminant Analysis,LDA)是一种监督学习的降维技术,与无监督的PCA不同的是,PCA是寻找数据集中方差最大的方向作为主成分分量的轴,而LDA是最优化分类的特征子空间。LDA的思想可以用一句话概括,就是“投影后类内方差最小,类间方差最大”。也就是说,要将数据在低维度上进行投影,投影后希望每一种类别数据的投影点尽可能的接近,而不同类别的数据的类别中心之间的距离尽可能的大。
       可能还是有点抽象,我们先看看最简单的情况。假设我们有两类数据 分别为红色和蓝色,如下图所示,这些数据特征是二维的,我们希望将这些数据投影到一维的一条直线,让每一种类别数据的投影点尽可能的接近,而红色和蓝色数据中心之间的距离尽可能的大。
在这里插入图片描述
       上图中提供了两种投影方式,哪一种能更好的满足我们的标准呢?从直观上可以看出,右图要比左图的投影效果好,因为右图的黑色数据和蓝色数据同类数据较为集中,且类别之间的距离明显,即满足类内方差最小,类间方差最大。左图则在边界处数据混杂。以上就是LDA的主要思想了,当然在实际应用中,我们的数据是多个类别的,我们的原始数据一般也是超过二维的,投影后的也一般不是直线,而是一个低维的超平面。

相关数学知识

1. 瑞利商(Rayleigh quotient)

       我们首先来看看瑞利商的定义。瑞利商是指这样的函数 R ( A , x ) R(A,x) R(A,x):
R ( A , x ) = x H A x x H x R(A, x) = \frac{x^HAx}{x^Hx} R(A,x)=xHxxHAx
       其中x为非零向量,而 A A A n × n n×n n×n H e r m i t a n Hermitan Hermitan矩阵。所谓的 H e r m i t a n Hermitan Hermitan矩阵就是满足共轭转置矩阵和自己相等的矩阵,即 A H = A AH=A AH=A。如果我们的矩阵A是实矩阵,则满足 A T = A A^T=A AT=A的矩阵即为 H e r m i t a n Hermitan Hermitan矩阵。
       具体的证明这里就不给出了。当向量 x x x是标准正交基时,即满足 x H x x^Hx xHx=1时,瑞利商退化为: R ( A , x ) = x H A x R(A,x)=x^HAx R(A,x)=xHAx,这个形式在谱聚类和PCA中都有出现。

2. 广义瑞利商(genralized Rayleigh quotient)

       广义瑞利商是指这样的函数 R ( A , B , x ) R(A,B,x) R(A,B,x):
R ( A , x ) = x H A x x H B x R(A,x)=\frac{x^HAx}{x^HBx} R(A,x)=xHBxxHAx
       其中 x x x为非零向量,而 A , B A,B A,B n × n n×n n×n H e r m i t a n Hermitan Hermitan矩阵。 B B B为正定矩阵。它的最大值和最小值是什么呢?其实我们只要通过将其通过标准化就可以转化为瑞利商的格式。我们令 x = B − 1 / 2 x ′ x=B^{−1/2}x^′ x=B1/2x,则分母转化为:
x H B x = x ′ H ( B − 1 / 2 ) H B B − 1 / 2 x ′ = x ′ H x ′ x^HBx=x^{′^H}(B^{−1/2})^HBB^{−1/2}x^′=x^{′^H}x^′ xHBx=xH(B1/2)HBB1/2x=xHx
       而分子转化为:
x H A x = x ′ H B − 1 / 2 A B − 1 / 2 x ′ x^HAx=x^{′^H}B^{−1/2}AB^{−1/2}x^′ xHAx=xHB1/2AB1/2x
       此时我们的 R ( A , B , x ) R(A,B,x) R(A,B,x)转化为 R ( A , B , x ′ ) R(A,B,x^′) R(A,B,x):
R ( A , B , x ′ ) = x ′ H B − 1 / 2 A B − 1 / 2 x ′ x ′ H x ′ R(A,B,x^′)=\frac{x^{′^H}B^{−1/2}AB^{−1/2}x^′}{x^{′^H}x^′} R(A,B,x)=xHxxHB1/2AB1/2x
       利用前面的瑞利商的性质,我们可以很快知道, R ( A , B , x ′ ) R(A,B,x^′) R(A,B,x)的最大值为矩阵 B − 1 / 2 A B − 1 / 2 B^{−1/2}AB^{−1/2} B1/2AB1/2的最大特征值,或者说矩阵 B − 1 A B^{−1}A B1A的最大特征值,而最小值为矩阵 B − 1 A B^{−1}A B1A的最小特征值。

LDA原理推导

        假设现在有数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x m , y m ) } D = \left \{(x_1, y_1), (x_2, y_2), ... , (x_m, y_m)\right \} D={(x1,y1),(x2,y2),...,(xm,ym)}, 其中任意样本 x i x_i xi n n n维向量, y i ∈ { C 1 , C 2 , . . . , C k } y_i \in \left \{C_1, C_2, ..., C_k\right \} yi{C1,C2,...,Ck}。我们定义 N j N_j Nj为第 j j j类样本的个数,其中 j ∈ { 1 , 2 , . . , k } j \in \left\{ 1, 2, .., k \right \} j{1,2,..,k} X j X_j Xj为第 j j j类样本的集合,而 μ j \mu_j μj为第 j j j类样本的均值向量,定义 Σ j ( j ∈ { 1 , 2 , . . . , k } ) \Sigma_j(j \in \left\{1, 2, ..., k\right\}) Σj(j{1,2,...,k})为第 j j j类样本的协方差矩阵。则:
μ j = 1 N j Σ x ∈ X j x \mu_j=\frac{1}{N_j} \Sigma_{x \in X_j}x μj=Nj1ΣxXjx
Σ j = Σ x ∈ X j ( x − μ j ) ( x − μ j ) T \Sigma_j=\Sigma_{x \in X_j}(x-\mu_j)(x-\mu_j)^T Σj=ΣxXj(xμj)(xμj)T
        假设我们的投影方向是向量 w w w,则对任意一个样本 x i x_i xi,它在方向w的投影是 w T x i w^T x_i wTxi,我们希望同一种类别的投影点尽可能接近,也就是类内协方差 ∑ i = 1 c ∑ j = 1 n i ∣ ∣ w T ( x j i − μ i ) ∣ ∣ 2 2 \sum_{i=1}^{c}\sum_{j=1}^{n_i}||w^T(x_j^i-\mu_i)||_2^2 i=1cj=1niwT(xjiμi)22 尽可能小。也就是说要最小化下面式子:
∑ i = 1 c ∑ j = 1 n i ∣ ∣ w T ( x j i − μ i ) ∣ ∣ 2 2 \sum_{i=1}^{c}\sum_{j=1}^{n_i}||w^T(x_j^i-\mu_i)||_2^2 i=1cj=1niwT(xjiμi)22
        另一方面,我们希望让不同类别的数据尽可能地远离,以便于我们分类。考虑类的中心,也就是均值 μ 1 , μ 2 , … , μ k μ_1,μ_2,…,μ_k μ1μ2,,μk只要最大化 ∑ i = 1 c n i ∣ ∣ w T ( μ i − μ ) ∣ ∣ 2 2 \sum_{i=1}^{c}n_i||w^T(\mu_i-\mu)||_2^2 i=1cniwT(μiμ)22
        综上所述,我们得到算法优化目标:
arg ⁡ max ⁡ J ( w ) = ∑ i = 1 c n i ∣ ∣ w T ( μ i − μ ) ∣ ∣ 2 2 ∑ i = 1 c ∑ j = 1 n i ∣ ∣ w T ( x j i − μ i ) ∣ ∣ 2 2 \arg\max J(w)=\frac{\sum_{i=1}^{c}n_i||w^T(\mu_i-\mu)||_2^2}{\sum_{i=1}^{c}\sum_{j=1}^{n_i}||w^T(x_j^i-\mu_i)||_2^2} argmaxJ(w)=i=1cj=1niwT(xjiμi)22i=1cniwT(μiμ)22
        为了表示方便,我们做如下定义:
S w = ∑ i = 1 c ∑ j = 1 n i ( x j i − μ i ) ( x j i − μ i ) T S_w=\sum_{i=1}^{c}\sum_{j=1}^{n_i}(x_j^i-\mu_i)(x_j^i - \mu_i)^T Sw=i=1cj=1ni(xjiμi)(xjiμi)T
S b = ∑ i = 1 c n i ( μ i − μ ) ( μ i − μ ) T S_b=\sum_{i=1}^{c}n_i(\mu_i-\mu)(\mu_i-\mu)^T Sb=i=1cni(μiμ)(μiμ)T
        优化目标重写为:
arg ⁡ max ⁡ J ( w ) = t r ( w T S b w ) t r ( w T S w w ) \arg\max J(w)=\frac{tr(w^TS_bw)}{tr(w^TS_ww)} argmaxJ(w)=tr(wTSww)tr(wTSbw)
        利用凸优化的知识,我们约束 w T S w w = 1 w^T S_w w=1 wTSww=1,这样使用拉格朗日乘子法就可以得到:
J ( w ) = w T S b w − λ ( w T S w w − 1 ) J(w)=w^TS_bw-\lambda (w^TS_ww-1) J(w)=wTSbwλ(wTSww1)
d J d w = 2 S b w − 2 λ S w w = 0 \frac{dJ}{dw}=2S_bw-2\lambda S_ww=0 dwdJ=2Sbw2λSww=0
        从而:
S b w = λ S w w S_bw=\lambda S_ww Sbw=λSww
        如果 S w S_w Sw可逆,就可以使用求解特征值分解的方法:
S w − 1 S b w = λ w S_w^{-1}S_bw=\lambda w Sw1Sbw=λw
        而在实际当中,S_w常常是不可逆的,我们可以通过以下两种方法解决:

  • S w = S w + γ I S_w=S_w + \gamma I Sw=Sw+γI,其中 γ \gamma γ是一个特别小的数,这样 S w S_w Sw一定可逆。
  • 先使用PCA对数据进行降维,使得在降维 后的数据上 S w S_w Sw可逆,再使用LDA求解。

matlab代码实现(以人脸识别为例)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%读取数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_face = load("ORL4646.mat").ORL4646;  % 所有人脸
face_train = [];  % 训练人脸
train_classes = []; % 训练人脸标记
face_test = [];  % 测试人脸
test_classes = []; % 测试人脸标记
% 一共40人,每个人取4张图片作为数据集
human_num = 40; % 数据集的人数
every_num = 10; % 每个人有几张照片
row = 46; % 照片行大小
col = 46; % 照片列大小
train_num = 5; % 用于训练的人脸数
test_num = 5; % 用于测试的人脸数
for i = 0 : human_num - 1
    for j = 1 : train_num
        tmp = all_face(:, :, i*every_num+j);
        new_face = reshape(tmp, row*col, 1); % 将像素矩阵按列拼接
        face_train = [face_train new_face]; % 加入到训练集中
    end
    for k = test_num : every_num
        tmp = all_face(:, :, i*every_num+k);
        new_face = reshape(tmp, row*col, 1); % 将像素矩阵按列拼接
        face_test = [face_test new_face]; % 加入到测试集中
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PCA降维%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PCAtrain_mean = mean(face_train,2); % PCA整体均值
PCAdimension = row*col;
PCA_mean = face_train(:, 1:train_num * human_num) - PCAtrain_mean;
[PCAV, PCAD] = eig(PCA_mean * PCA_mean');
[PCAd, PCAind] = sort(diag(PCAD), "descend");
Ds = PCAD(PCAind, PCAind);
Vs = PCAV(:, PCAind);
A = Vs(:, 1:PCAdimension);
LDAface_train =  face_train;
LDAface_test = face_test;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算Sb和Sw%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LDAtrain_mean = mean(LDAface_train,2); % LDA整体均值
Sb = 0; % 类间散度矩阵
Sw = 0; % 类内散度矩阵
N_j = train_num; % 每个类的样本数
avcolumn = zeros(PCAdimension,PCAdimension);  % 计算每个类的均值                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
for i = 1:human_num
    avcolumn(:,i) = mean(LDAface_train(:,(i-1)*N_j+1:i*N_j),2);
end
for i = 0 : (human_num - 1)
    for k = 1 : N_j
        temp2 = LDAface_train(:, i*N_j+k) - avcolumn(:,i+1); % 用于计算类内散度
        Sw = Sw + temp2 * temp2';
    end
end
St = zeros(PCAdimension,PCAdimension); % 计算St
for i = 1:N_j*human_num
    St = St + (LDAface_train(:,i)-LDAtrain_mean)*(LDAface_train(:,i)-LDAtrain_mean)';
end
Sb = St - Sw;
%[U, Sigma, DD] = svd(Sw);
%Sigma_diag = diag(Sigma);
%Sigma = Sigma + (1e-11)*eye(size(Sigma, 1), size(Sigma, 2));
%Sw = U * Sigma * DD;
Sw = Sw + (1e-7)*eye(size(Sw, 1), size(Sw, 2));
%%%%%%%%%%%%%%%%方法一%%%%%%%%%%%%%%%%%%%%
[V, D] = eig(Sw\Sb);% 将特征向量按照特征值降序
[LDAd, LDAind] = sort(diag(D), "descend");

%%%%%%%%%%%%%%%%方法二%%%%%%%%%%%%%%%%%%%%
%Sb = Sb + (1e-7)*eye(size(Sb, 1), size(Sb, 2));
%[V, D] = eig(Sb\Sw);% 将特征向量按照特征值降序
%[LDAd, LDAind] = sort(diag(D));

%%%%%%%%%%%%%%%%方法三%%%%%%%%%%%%%%%%%%%%
%[V, D] = eig(Sb - Sw);% 将特征向量按照特征值降序
%[LDAd, LDAind] = sort(diag(D),"descend");

Ds = D(LDAind, LDAind);
LDAeigenVector = V(:, LDAind);
%dimension = 39;
%W = eigenVector(:, 1 :dimension); % 获取投影矩阵的转置
%FisherFaces = W' * face_train;
dimensions = [];
rates = [];
for d = 3:3
    dimension = d * 1;  % 降维后的维数
    dimensions = [dimensions dimension]; 
    W = LDAeigenVector(:, 1 :dimension); 
    FisherFaces = W' * LDAface_train;
    TestFaces = W' * LDAface_test;
    for i = 0 : 2
        tmp = FisherFaces(:, (i*10+1):(i+1)*10);
        %tmp = FisherFaces(:, (i*10+1):(i*10+1));
        color = rand(1, 3);
        plot3(tmp(1,:),tmp(2,:),tmp(3,:),'marker','.','linestyle','none','color',color, 'MarkerSize',10);
        hold on;
    end
    hold on;
    axes("Position",[0.58 0.45 0.05 0.05]);imshow(reshape(face_train(:,1),46,46),[]);
    hold on;
    axes("Position",[0.3 0.68 0.05 0.05]);imshow(reshape(face_train(:,11),46,46),[]);
    hold on;
    axes("Position",[0.3 0.48 0.05 0.05]);imshow(reshape(face_train(:,21),46,46),[]);
    %subplot(7, 7, d)
    %imshow(reshape(W*FisherFaces(:, 1), row, col), []);
    correct_sum = 0;
    for test_index = 1 : size(TestFaces, 2)
        test_face = TestFaces(:, test_index);
        Euc_dist = [];
        for i = 1 : size(FisherFaces, 2)
            tmp = ( norm(test_face - FisherFaces(:, i), 2))^2; % 计算L2范数
            Euc_dist = [Euc_dist tmp];  % 用于存放距离
        end
        [Euc_dist_min , Recognized_index] = min(Euc_dist); % 获得最小距离
        disp([Recognized_index test_index]);
        if ceil(test_index/test_num) == ceil(Recognized_index/train_num)  % 进行分类
            correct_sum = correct_sum + 1;
        end
    end
    rate = correct_sum / size(TestFaces, 2); % 计算准确率
    rates = [rates rate];
end
% 画折线图
%disp(rates);
%p = plot(dimensions, rates,'k-o','linewidth', 2, 'markersize', 4);
%p.Color = "blue";
%grid on;
%xlabel('维度','fontname','标楷体', 'fontweight', 'bold', 'fontsize', 12 );
%ylabel('识别率', 'fontname', '标楷体', 'fontweight', 'bold', 'fontsize',12);
%xlim([1, human_num-1]);
%ylim([0, 1]);

Python代码实现(以人脸识别为例)

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsClassifier
from sklearn import metrics
from mpl_toolkits.mplot3d import Axes3D

def getdata(eachNum):
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\ORL4646.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['ORL4646']  # 删除作者版本等信息,只采用人脸数据
    datamatrix = np.reshape(datamatrix, (46 * 46, 400))  # 将人脸数据转化为二维数组
    datamatrix = np.array(datamatrix)
    datamatrix = np.transpose(datamatrix) #
    train_face = [];test_face = []; train_label = [];test_label = []
    for i in range(40):
        for j in range(10):
            if(j<eachNum):
                train_face.append(datamatrix[i*10+j])
                train_label.append(i+1)
            else:
                test_face.append(datamatrix[i*10+j])
                test_label.append(i+1)
    train_face = np.asarray(train_face)
    train_label = np.asarray(train_label)
    test_face = np.asarray(test_face)
    test_label = np.asarray(test_label)
    '''
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\YaleB_32x32.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['fea']
    train_face = [];test_face = [];train_label = [];test_label = []
    for i in range(38):
        for j in range(10):
            if (j < eachNum):
                train_face.append(datamatrix[i * 64 + j])
                train_label.append(i + 1)
            else:
                test_face.append(datamatrix[i * 64 + j])
                test_label.append(i + 1)
    train_face = np.asarray(train_face)
    train_label = np.asarray(train_label)
    test_face = np.asarray(test_face)
    test_label = np.asarray(test_label)
    '''
    '''
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\Yale5040165.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['Yale5040165']
    datamatrix = np.reshape(datamatrix, (50 * 40, 165))  # 将人脸数据转化为二维数组
    datamatrix = np.array(datamatrix)
    datamatrix = np.transpose(datamatrix)
    train_face = [];test_face = [];train_label = [];test_label = []
    for i in range(15):
        for j in range(11):
            if (j < eachNum):
                train_face.append(datamatrix[i * 11 + j])
                train_label.append(i + 1)
            else:
                test_face.append(datamatrix[i * 11 + j])
                test_label.append(i + 1)
    train_face = np.asarray(train_face)
    train_label = np.asarray(train_label)
    test_face = np.asarray(test_face)
    test_label = np.asarray(test_label)
    return train_face,train_label,test_face,test_label,datamatrix
    '''
    return train_face,train_label,test_face,test_label,datamatrix

def getdataByRecog(eachNum):
    '''
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\ORL4646.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['ORL4646']  # 删除作者版本等信息,只采用人脸数据
    datamatrix = np.reshape(datamatrix, (46 * 46, 400)) #将人脸数据转化为二维数组
    datamatrix = np.array(datamatrix)
    datamatrix = np.transpose(datamatrix)
    trainMatrix = []
    testMatrix = []
    for i in range(40): #ORL共40个人
        for j in range(10): #每个人10张照片
            if j < eachNum: #每个人用n张照片进行训练
                trainMatrix.append(datamatrix[i*10 + j])
            else: #剩下的进行测试
                testMatrix.append(datamatrix[i*10 + j])
    trainMatrix = np.array(trainMatrix)
    testMatrix = np.array(testMatrix)
    #print(trainMatrix.shape)
    #print(testMatrix.shape)
    return trainMatrix,testMatrix
    :param eachNum:
    :return:
    '''
    '''
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\YaleB_32x32.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['fea']
    trainMatrix = []
    testMatrix = []
    for i in range(38):  # YaleB共个人
        for j in range(10):  # 每个人64张照片
            if j < eachNum:  # 每个人用n张照片进行训练
                trainMatrix.append(datamatrix[i * 64 + j])
            else:  # 剩下的进行测试
                testMatrix.append(datamatrix[i * 64 + j])
    trainMatrix = np.array(trainMatrix)
    testMatrix = np.array(testMatrix)
    return trainMatrix, testMatrix
    '''
    dataFile = 'D:\大学作业\大二下机器学习\算法常用数据\Yale5040165.mat'
    data = scio.loadmat(dataFile)
    datamatrix = data['Yale5040165']
    datamatrix = np.reshape(datamatrix, (50 * 40, 165))  # 将人脸数据转化为二维数组
    datamatrix = np.array(datamatrix)
    datamatrix = np.transpose(datamatrix)
    trainMatrix = []
    testMatrix = []
    for i in range(15):  # Yale共15个人
        for j in range(11):  # 每个人张照片
            if j < eachNum:  # 每个人用n张照片进行训练
                trainMatrix.append(datamatrix[i * 11 + j])
            else:  # 剩下的进行测试
                testMatrix.append(datamatrix[i * 11 + j])
    trainMatrix = np.array(trainMatrix)
    testMatrix = np.array(testMatrix)
    return trainMatrix, testMatrix

def LDA(face,label,k):
    x,y = face.shape
    classes = np.unique(label)
    meanAll = face.mean(axis=0)
    Sb = np.zeros((y, y),dtype = np.float32)
    Sw = np.zeros((y, y), dtype=np.float32)
    for i in classes:
        face_i = face[np.where(label == i)[0],:]
        mean_i = face_i.mean(axis = 0)
        n = face_i.shape[0]
        Sw = Sw + np.dot((face_i - mean_i).T,(face_i-mean_i))
        Sb = Sb + n*np.dot((mean_i - meanAll).T,(mean_i - meanAll))
    #print("Sw");print(Sw);print("end");print("Sb");print(Sb);print("end")
    tmp = np.ones(Sw.shape)/500
    Sw = tmp + Sw
    matrix = np.linalg.inv(Sw) @ Sb
    eigenvalue,eigenvector = np.linalg.eigh(matrix)
    #print("vec");print(eigenvector);print("end")
    index = np.argsort(eigenvalue)  # 排序
    index = index[:-(k + 1):-1]
    select = eigenvector[:, index]
    return select

def pca(matrix,k):
    matrix = np.float32(np.mat(matrix))
    meanMatrix = np.mean(matrix,axis=0)
    matrix = matrix - meanMatrix #去平均值
    covMatrix = np.cov(matrix.T,rowvar=1) #协方差矩阵
    eigenvalue,eigenvector = np.linalg.eigh(np.mat(covMatrix)) #特征值和特征向量
    index = np.argsort(eigenvalue) #排序
    index = index[:-(k+1):-1]
    select = eigenvector[:,index] #最大的特征值对应的特征向量
    return select.T

def pca_recog(matrix,k):
    matrix = np.float32(np.mat(matrix))
    eachMean = []
    #eachMean = np.array(eachMean)
    for i in range(matrix.shape[0]):
        eachMean.append(matrix[i].mean(axis=0))
    eachMean = np.array(eachMean)
    eachMean = eachMean.reshape(eachMean.shape[0],eachMean.shape[2])
    matrix = matrix.T
    meanMatrix = matrix.mean(axis=1)  # 按行计算平均
    matrix = matrix - meanMatrix  # 去平均值
    u,sigma,v = np.linalg.svd(matrix)
    w = u[:,:k]
    w = np.array(w)
    eachMean = np.array(eachMean)
    return w,eachMean

def recognize_pca():
    eachNum = 5 #每个人用来训练的照片数
    trainMatrix,testMatrix = getdataByRecog(eachNum)
    #print(testMatrix.shape)
    trainNum = trainMatrix.shape[0]
    testNum = testMatrix.shape[0]
    acc = []
    dimen = []
    w,eachMean = pca_recog(trainMatrix,180)
    project = []
    for each in eachMean:
        tmp = np.dot(w.T,each.T)
        project.append(tmp.T)
    project = np.array(project)
    correct = 0
    for i in range(testNum):
        tmp = np.reshape(testMatrix[i],(50*40,1))
        test = np.dot(w.T , tmp)
        disMatrix = []
        for j in range(project.shape[0]):
            distance = np.linalg.norm(test.T - project[j])
            disMatrix.append(distance)
        classes = int(disMatrix.index(min(disMatrix))/eachNum) #eachNum张为一类
        if classes == int(i/(10-eachNum)): #10-eachNum张为一类
            correct += 1
    accuracy = float(correct) / testNum
    print(accuracy)

def KNN_classify(face,train_face,train_label,k):
    dis = np.sum(np.power((train_face - face),2),axis=1)
    index = np.argsort(dis)[:k]
    weight = []
    for i in range(k):
        weight.append(dis[index[k-1]] - dis[index[i]] / (dis[index[k-1]] - dis[index[0]]))
    count = [0 for i in range (train_face.shape[0])]
    tmp = 0
    for i in index:
        count[train_label[i]] += 1 + weight[tmp]
        tmp += 1
    label = np.argmax(count);
    return label

def KNN(train_faces, train_labels,test_faces, test_labels,k):
    sum = test_faces.shape[0]
    err = 0
    for i in range(sum):
        count = KNN_classify(test_faces[i],train_faces,train_labels,k)
        if count != test_labels[i]:
            err += 1
    acc = (sum - err) / sum
    return acc

def view():
    eachNum = 7  # 每个人拿出5张照片进行训练
    train_face, train_label, test_face, test_label, face = getdata(eachNum)

    W = LDA(train_face, train_label, 3, eachNum)
    # print(W)
    project = np.dot(W, face.T).T
    print(project.shape)
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    colors = ['b', 'g', 'r']
    for i in range(3):
        tmpx = []
        tmpy = []
        tmpz = []
        for j in range(10):
            tmpx.append(project[10 * i + j][0])
            tmpy.append(project[10 * i + j][1])
            tmpz.append(project[10 * i + j][2])
        ax.scatter(tmpx, tmpy, tmpz, c=colors[i])
    # ax.title("LDA")
    plt.show()

    matrixW = pca(train_face, 2)
    matrixW = np.array(matrixW)
    # print(matrixW)
    project_pca = np.dot(matrixW, face.T).T
    # print(project_pca.shape)
    fig = plt.figure()
    ax1 = fig.add_subplot(111, projection='3d')
    colors = ['b', 'g', 'r']
    for i in range(3):
        tmpx = []
        tmpy = []
        tmpz = []
        for j in range(10):
            tmpx.append(project_pca[10 * i + j][0])
            tmpy.append(project_pca[10 * i + j][1])
            tmpz.append(project_pca[10 * i + j][1])
        ax1.scatter(tmpx, tmpy, tmpz, c=colors[i])
    # ax1.title("PCA")
    plt.show()

def rebuild():
    eachNum = 7  # 每个人拿出5张照片进行训练
    train_face, train_label, test_face, test_label, face = getdata(eachNum)
    print(face[1].shape)
    W = LDA(train_face, train_label, 39, eachNum)
    lowmatrix = np.mat(face[10]) * W.T # 1,2116 * 2116*d
    fisherface = lowmatrix*W #1,d * d*2116
    fisherface = np.reshape(fisherface,(46,46))
    plt.figure()
    plt.axis('off')
    plt.imshow(fisherface, cmap=plt.cm.gray)
    plt.show()

    W_pca = pca(train_face, 180)
    lowmatrix_pca = np.mat(face[10]) * W_pca.T
    eigenface = lowmatrix_pca * W_pca
    eigenface = np.reshape(eigenface,(46,46))
    plt.figure()
    plt.axis('off')
    plt.imshow(eigenface, cmap=plt.cm.gray)
    plt.show()

    plt.figure()
    plt.axis('off')
    plt.imshow(face[10].reshape(46,46), cmap=plt.cm.gray)
    plt.show()

def recognize_lda():
    eachNum = 3  # 每个人拿出eachNum张照片进行训练
    train_face, train_label, test_face, test_label, face = getdata(eachNum)

    W = LDA(train_face, train_label, 39)
    train = np.dot(train_face,W)
    test = np.dot(test_face,W)
    acc = KNN(train,train_label,test,test_label,3)
    print("acc=",acc);

if __name__ == "__main__":
    #view()
    #rebuild()
    #recognize_lda()
    recognize_pca()

相关证明: S t = S w + S b S_t=S_w+S_b St=Sw+Sb

在这里插入图片描述

LDA主要优点

  1. 在降维过程中可以使用类别的先验知识经验,而像PCA这样的无监督学习则无法使用类别先验知识。
  2. LDA在样本分类信息依赖均值而不是方差的时候,比PCA之类的算法较优。

LDA主要缺点

  1. LDA不适合对非高斯分布样本进行降维,PCA也有这个问题。
  2. LDA降维最多降到类别数k-1的维数,如果我们降维的维度大于k-1,则不能使用LDA。当然目前有一些LDA的进化版算法可以绕过这个问题。
  3. LDA在样本分类信息依赖方差而不是均值的时候,降维效果不好。
  4. LDA可能过度拟合数据。

写在最后

        本文主要为个人学习的成果总结,部分内容参考网上相关资料,如有侵权联系删除。另外,也欢迎各位对本文进行转载,转载请附上原文链接。读者在阅读本文过程中如有疑问,欢迎随时交流。

现在我们回到LDA原理上,我们在第一节说讲到了LDA希望投影后希望同一种类别数据的投影点尽可能的接近,而不同类别的数据的类别中心之间的距离尽可能的大,但是这只是一个感官的度量。现在我们首先从比较简单的二类LDA入手,严谨的分析LDA原理。     假设我们的数据集D={(x1,y1),(x2,y2),...,((xm,ym))}D={(x1,y1),(x2,y2),...,((xm,ym))},其中任意样本xixi为n维向量,yi∈{0,1}yi∈{0,1}。我们定义Nj(j=0,1)Nj(j=0,1)为第j类样本的个数,Xj(j=0,1)Xj(j=0,1)为第j类样本的集合,而μj(j=0,1)μj(j=0,1)为第j类样本的均值向量,定义Σj(j=0,1)Σj(j=0,1)为第j类样本的协方差矩阵(严格说是缺少分母部分的协方差矩阵)。     μjμj的表达式为: μj=1Nj∑x∈Xjx(j=0,1) μj=1Nj∑x∈Xjx(j=0,1)     ΣjΣj的表达式为: Σj=∑x∈Xj(x−μj)(x−μj)T(j=0,1) Σj=∑x∈Xj(x−μj)(x−μj)T(j=0,1)     由于是两类数据,因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量ww,则对任意一个样本本xixi,它在直线ww的投影为wTxiwTxi,对于我们的两个类别的中心点μ0,μ1μ0,μ1,在在直线ww的投影为wTμ0wTμ0和wTμ1wTμ1。由于LDA需要让不同类别的数据的类别中心之间的距离尽可能的大,也就是我们要最大化||wTμ0−wTμ1||22||wTμ0−wTμ1||22,同时我们希望同一种类别数据的投影点尽可能的接近,也就是要同类样本投影点的协方差wTΣ0wwTΣ0w和wTΣ1wwTΣ1w尽可能的小,即最小化wTΣ0w+wTΣ1wwTΣ0w+wTΣ1w。综上所述,我们的优化目标为: argmaxwJ(w)=||wTμ0−wTμ1||22wTΣ0w+wTΣ1w=wT(μ0−μ1)(μ0−μ1)TwwT(Σ0+Σ1)w argmax⏟wJ(w)=||wTμ0−wTμ1||22wTΣ0w+wTΣ1w=wT(μ0−μ1)(μ0−μ1)TwwT(Σ0+Σ1)w     我们一般定义类内散度矩阵SwSw为: Sw=Σ0+Σ1=∑x∈X0(x−μ0)(x−μ0)T+∑x∈X1(x−μ1)(x−μ1)T Sw=Σ0+Σ1=∑x∈X0(x−μ0)(x−μ0)T+∑x∈X1(x−μ1)(x−μ1)T     同时定义类间散度矩阵SbSb为: Sb=(μ0−μ1)(μ0−μ1)T Sb=(μ0−μ1)(μ0−μ1)T     这样我们的优化目标重写为: argmaxwJ(w)=wTSbwwTSww argmax⏟wJ(w)=wTSbwwTSww     仔细一看上式,这不就是我们的广义瑞利商嘛!这就简单了,利用我们第二节讲到的广义瑞利商的性质,我们知道我们的J(w)J(w)最大值为矩阵S−12wSbS−12wSw−12SbSw−12的最大特征值,而对应的ww为S−12wSbS−12wSw−12SbSw−12的最大特征值对应的特征向量! 而S−1wSbSw−1Sb的特征值和S−12wSbS−12wSw−12SbSw−12的特征值相同,S−1wSbSw−1Sb的特征向量w′w′和S−12wSbS−12wSw−12SbSw−12的特征向量ww满足w′=S−12www′=Sw−12w的关系!     注意到对于二类的时候,SbwSbw的方向恒为μ0−μ1μ0−μ1,不妨令Sbw=λ(μ0−μ1)Sbw=λ(μ0−μ1),将其带入:(S−1wSb)w=λw(Sw−1Sb)w=λw,可以得到w=S−1w(μ0−μ1)w=Sw−1(μ0−μ1), 也就是说我们只要求出原始二类样本的均值和方差就可以确定最佳的投影方向ww了。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Kim‘s blog

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值