以下使用数据集为:COIL20.mat(1440X1024)

下载链接:https://pan.baidu.com/s/13kXvS4elgHmuuttS8arFqA

提取码:rvyo

LDA介绍:

LDA是一种常用的数据降维方式，它属于监督式降维（特征提取）。通常在分类之前降低数据维度，使得不同类之间的数据更加的泾渭分明。

原理解释：

基本原理：将高维的样本数据投影到最佳判别向量空间,以达到特征提取(维数约简)的效果,投影后保证样本数据在新的子空间有最大的类间距离和最小的类内距离,即在该子空间中具有最佳的可分离性。

数学解释: 通过一个线性变换,将(为的矩阵)中的样本数据映射到(从d维降到k维),且希望该变换将属于同一类的样本映射得越近越好(即最小的类内距离)，而将不同类的样本映射得越远越好 (即最大的类间距离)。同时还能尽能多地保留样本数据的判别信息。

我们所要做的就是找到一个线性变换,满足上述描述。

更加精确的描述就是： 显然 w应该是一个 的矩阵。保证矩阵

优化目标：

有了问题，我们就要知道如何让我们找到的答案越来越好。

先从简单的二分类问题出发：

我们尝试将二维点坐标映射到一条直线上，显然我们要选择右图中的映射方式：

• 映射后，不同类之间有着明显的分界
• 同类型数据分布更加紧凑

我们尝试找出两个数学表达表达上述两种条件。

设置两类样本集合,分别对应样本均值   

经过线性变化后，得到映射到直线上的。()

这样我们可以将两类中心之间的距离近似看为类间距离。

• 之间距离越大越好（类间离散度）

同样我们要求类内距离小(同类数据方差小)

• 内数据集中在附近(类内离散度)

           ，其中

所以最后的优化目标函数是：

更进一步的化简：

注：多分类下，，是总体均值

最终转化为:

...

更详细的化简可以翻阅《模式识别》一书中的LDA线性判别分析。

最终化简为：

(这是在可逆的情况下)。

（如果矩阵不可逆，可为奇异矩阵加上一个单位矩阵，单位矩阵乘以一个很小的系数，如0.001)

结论:最优w就是的特征向量矩阵，而这个公式也称为Fisher线性判别

注意：如果我们要求d维降到k维，则我们的w是最大的k个特征向量构成的矩阵。

也许，你求出的w矩阵是复数矩阵，那么取w得实部。

代码实现：

# -*- coding: utf-8 -*-
"""
Created on Wed May 15 14:48:32 2019
Python 3.6.8
@author: yh
"""
from numpy import linalg as LA
from scipy import io
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
import tensorflow as tf
"""
LDA线性判别分析
数据降维+KNN分类

"""

np.random.seed(0) 
"""
导入数据
这个数据集1400X1024
20X72
20中类别,每类72个32x32的灰度图像
"""
data=io.loadmat('COIL20.mat')
x=data['fea']
y=data['gnd']


def Std(Sw):
    u=np.mean(Sw)
    s=np.std(Sw)
    Sw=(Sw-u)/s
    return Sw

#X是数据集,k是降到k维,n是数据集中类的个数,m是每一个样本的维度
def Fisher(X,k,n,m):
    #求每一类的均值
    mu=np.zeros((n,1,m))
    for i in range(n):
        for j in range(len(X[i])):
            mu[i]=mu[i]+X[i][j]
        mu[i]=mu[i]/len(X[i])
    #总体均值
    u=np.mean(mu,axis=0)
    Sw=np.zeros((m,m))
    Sb=np.zeros((m,m))
    for i in range(n):
        for j in range(len(X[i])):
            #注意np.dot才是矩阵乘法,‘*’对于array来说是点乘
            Sw=Sw+np.dot((X[i][j]-mu[i]).T,(X[i][j]-mu[i]))
    for i in range(n):
        Sb=Sb+np.dot((mu[i]-u).T,(mu[i]-u))
    #单位矩阵
    I=np.identity(m)
    #归一化
    Sw=Std(Sw);
    Sw=Sw+0.000000001*I
    val,vec=np.linalg.eig(np.dot(np.mat(Sw).I ,Sb))
    #求最大K个特征向量编号
    index=np.argsort(-val)
    #最大K个特征向量所组成的矩阵
    w=vec[index[0:k]]
    #取w矩阵实部返回
    return w.real

    

def Class_KNN(X,Y):
    #KNN分类
    #自动分出训练集,测试集
    x_train,x_test,y_train,y_test = train_test_split(X,Y)
    knn = KNeighborsClassifier() 
    knn.fit(x_train,y_train.ravel())
    accuracy=knn.score(x_test,y_test,sample_weight=None)
    return accuracy
    

def Showimage(X,x,n,m):
    ###这里可以修改最大降维维度
    k=100
    K=np.arange(1,k)
    Ac=np.zeros((1,k-1))
    print(K.shape,Ac.shape)
    for i in range(len(K)):
        w=Fisher(X,K[i],n,m)
        z=np.dot(x,w.T)
        Ac[0][i]=Class_KNN(z,y)
    #画出k-ac图片
    plt.plot(K,Ac.ravel(),'-r')
    plt.xlabel('K-Dimesions')
    plt.ylabel('Accuracy')


if __name__=="__main__":
        n,m=x.shape
        print(n,m)
        X=np.zeros((n//72,72,m))
        for i in range(n//72):
            xi=x[i*72:i*72+72,:]
            X[i]=xi
        Showimage(X,x,n//72,m)
    
            

结果展示: