降维：LDA推导及iris实例

线性判别分析

\qquad 线性判别分析在降维中的运用：用过将高维的样本投影到最佳的矢量空间，投影后的样本在新的子空间有最大类间距离以及最小类内距离。也就是使投影后的数据内间散度矩阵大，类内散度矩阵小，等价于相同类别数据紧凑，不用类别离得远。

PCA与LDA之间的区别？

PCA	LDA
无监督的数据降维方法	有监督的数据降维方法
不使用数据标签	使用用数据标签
投影到方差最大的相互正交的方向上（保留更多的信息）	尽可能多的保留样本信息；寻找使样本尽可能好分的方向；同类近，不同尽可能远
	最多降到类别数N-1维
	可用于分类

LDA推导

二分类

数据集：D={($x_1$,$y_1$),($x_2$,$y_2$)…($x_m$,$y_m$)} \qquad $y_i$∈{0,1} \qquad $x_i$∈$R^{n\times 1}$
n：特征数
m：样本数量
$m_i$：第i类样本数量
待求：投影方向$w$ \qquad $w$∈$R^{n\times 1}$
${\mu }_i$：第i类数据集的样本均值
${\sum }_i$:第类数据集变换前的散度矩阵
\qquad \qquad ${\sum }_i$=${\sum }_{x\in X_i}$($x$-${\mu }_i$)($x$-${\mu }_i$${)}^T$
 
$X_i$：第i类拥有的数据集
$w^T$${\mu }_i$：第i类中心在$w$上的投影 \qquad $w^T$${\mu }_i$∈$R^{1\times 1}$
$w^T$$x_i$：第i个样本在$w$上的投影 \qquad $w^T$$x_i$∈$R^{1\times 1}$
$w^T$${\sum }_i$ $w$：投影后的方差 \qquad $w^T$${\sum }_i$ $w$∈$R^{1\times 1}$
\qquad \qquad ${\sum }_{x\in X_i}$($w^T$$x$-$w^T$${\mu }_i$)($w^T$$x$-$w^T$${\mu }_i$${)}^T$
\qquad \qquad =${\sum }_{x\in X_i}$$w^T$($x$-${\mu }_i$)($x$-${\mu }_i$${)}^T$$w$
\qquad \qquad =$w^T$${\sum }_{x\in X_i}$($x$-${\mu }_i$)($x$-${\mu }_i$${)}^T$$w$
\qquad \qquad =$w^T$${\sum }_i$$w$
 
类内方差小，类间距离大
 
目标函数：
\qquad \qquad \qquad max J($w$)=$\frac{||w^T{\mu }_0-w^T{\mu }_1|{|}^2}{w^T{\sum }_{0}w+w^T{\sum }_{1}w}$
 
\qquad \qquad \qquad \qquad              \;\;\;\;\;\; =$\frac{||w^T{\mu }_0-w^T{\mu }_1|{|}^2}{w^T{\sum }_{0}w+w^T{\sum }_{1}w}$
 
\qquad \qquad \qquad \qquad              \;\;\;\;\;\; =$\frac{||w^T({\mu }_0-{\mu }_1)|{|}^2}{w^T({\sum }_0+{\sum }_1)w}$
 
\qquad \qquad \qquad \qquad              \;\;\;\;\;\; =$\frac{w^T({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w}{w^T({\sum }_0+{\sum }_1)w}$
 
类内散度矩阵（半正定矩阵，对称阵）：$S_w$=${\sum }_0$+${\sum }_1$
类间散度矩阵（对称阵）：$S_b$=(${\mu }_0$-${\mu }_1$)(${\mu }_0$-${\mu }_1$${)}^T$
 
目标函数转化为：
\qquad \qquad \qquad max J($w$)=$\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$
 
由广义瑞丽商得:J($w$)的最大值为$S_w^{-1}$$S_b$的最大特征值
 
求投影方向$w$:最重要的是确定$w$的方向
此时目标函数优化为：max $w^T$$S_b$$w$              \;\;\;\;\;\; s.t.$w^T$$S_w$w$=1
拉格朗日：
\qquad \qquad \qquad L($w$，λ)=$w^T$$S_b$$w$-λ($w^T$$S_w$w$-1)L对$w$求导，并且让其等于0
\qquad \qquad \qquad （$S_b$+$S_b^T$）$w$-λ（$S_w$+$S_w^T$）$w$=0
\qquad \qquad \qquad \qquad \qquad $S_b$$w$=λ$S_w$$w$ \qquad \qquad \qquad （1）
又
\qquad $S_b$$w$=(${\mu }_0$-${\mu }_1$)(${\mu }_0$-${\mu }_1$${)}^T$$w$

且$S_b$$w$的方向恒为${\mu }_0$-${\mu }_1$

假设 \qquad $S_b$$w$=λ(${\mu }_0$-${\mu }_1$) \qquad \qquad \qquad （2）

即让 \qquad λ=(${\mu }_0$-${\mu }_1$${)}^T$$w$

（2）带入（1）得：
\qquad \qquad \qquad $w$=$S_w^{-1}$(${\mu }_0$-${\mu }_1$)

多分类

数据集：D={($x_1$,$y_1$),($x_2$,$y_2$)…($x_m$,$y_m$)} \qquad $y_i$∈{1,2…N} \qquad $x_i$∈$R^{n\times 1}$
N：类数
n：特征数
m：样本数量
$m_i$：第i类样本数量
待求：投影矩阵$W$={$w_1$,$w_2$…$w_d$} \qquad $w$∈$R^{n\times d}$
${\mu }_i$：第i类数据集的样本均值
$\mu$：整个数据集的样本均值
 
全局散度矩阵：$S_t$=$S_b$+$S_w$
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$${\sum }_{x\in X_i}$($x$-$\mu$)($x$-$\mu$${)}^T$

类内散度矩阵：$S_w$=${\sum }_{i=1}^N$${\sum }_{x\in X_i}$($x$-${\mu }_i$)($x$-${\mu }_i$${)}^T$

类间散度矩阵：$S_b$=$S_t$-$S_w$
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$${\sum }_{x\in X_i}$[($x$-$\mu$)($x$-$\mu$${)}^T$-($x$-${\mu }_i$)($x$-${\mu }_i$${)}^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$${\sum }_{x\in X_i}$[$x$$x^T$-$x$${\mu }^T$-$\mu$$x^T$+$\mu$${\mu }^T$-$x$$x^T$+$x$${\mu }_i^T$+${\mu }_i$$x^T$-${\mu }_i$${\mu }_i^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$${\sum }_{x\in X_i}$[-$x$${\mu }^T$-$\mu$$x^T$+$\mu$${\mu }^T$+$x$${\mu }_i^T$+${\mu }_i$$x^T$-${\mu }_i$${\mu }_i^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$[-$m_i$${\mu }_i$${\mu }^T$-$m_i$$\mu$${\mu }_i^T$+$m_i$$\mu$${\mu }^T$+$m_i$${\mu }_i$${\mu }_i^T$+$m_i$${\mu }_i$${\mu }_i^T$-$m_i$${\mu }_i$${\mu }_i^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$[-$m_i$${\mu }_i$${\mu }^T$-$m_i$$\mu$${\mu }_i^T$+$m_i$$\mu$${\mu }^T$+$m_i$${\mu }_i$${\mu }_i^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$$m_i$[-${\mu }_i$${\mu }^T$-$\mu$${\mu }_i^T$+$\mu$${\mu }^T$+${\mu }_i$${\mu }_i^T$]
\qquad \qquad \qquad        \;\;\; =${\sum }_{i=1}^N$$m_i$(${\mu }_i$-$\mu$)(${\mu }_i$-$\mu$${)}^T$
 
目标函数：
\qquad \qquad argmax J($w_i$)=$\frac{{\sum }_{i=1}^d{w}_i^T{S}_b{w}_i}{{\sum }_{i=1}^d{w}_i^T{S}_w{w}_i}$
\qquad \qquad \qquad \qquad          \;\;\;\; =${\sum }_{i=1}^d$$\frac{w_i^T{S}_b{w}_i}{w_i^T{S}_w{w}_i}$
由广义瑞丽熵得：J($w_i$)的最大值为$S_w^{-1}$$S_b$的最大d个特征值之和

求$w_i$：
\qquad \qquad argmax J($w_i$)=${\sum }_{i=1}^d$$\frac{w_i^T{S}_b{w}_i}{w_i^T{S}_w{w}_i}$
\qquad \qquad \qquad \qquad          \;\;\;\; =$\frac{tr(W^T{S}_{b}W)}{tr(W^T{S}_{w}W)}$

等价于：
\qquad \qquad argmax tr($W^T$$S_b$$W$) \qquad \qquad s.t. tr($W^T$$S_w$$W$)=1
拉格朗日：
\qquad L（W,∧）=-tr($W^T$$S_b$$W$)+tr（∧（$W^T$$S_w$$W$-I））
∧=diag(${\lambda }_1$,${\lambda }_2$…${\lambda }_d$)
L对W求导，并让其为0：
\qquad \qquad \qquad （$S_b$+$S_b^T$）$W$-∧（$S_w$+$S_w^T$）$W$=0
\qquad \qquad \qquad \qquad \qquad $S_b$$W$=∧$S_w$$W$
对任意${\lambda }_i$都有：
\qquad \qquad \qquad \qquad \qquad $S_b$$W$=${\lambda }_i$$S_w$$W$
=> \qquad \qquad \qquad \qquad          \;\;\;\; $S_w^{-1}$$S_b$$W$=${\lambda }_i$$W$

$w_i$为$S_w^{-1}$$S_b$的d个最大特征值对应的特征向量

iris实例

LDA降维步骤

1. 计算每类的均值向量
2. 计算整个样本的均值向量
3. 计算类内散度矩阵$S_w$以及类间散度矩阵$S_b$
4. 计算$S_w^{-1}$$S_b$的特征值以及特征向量，选取前d个最大特征值对应的特征向量，组成矩阵W
5. 对原始数据降维。降维数据Y=$W^T$X

单步

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

iris=datasets.load_iris()
#一共三类数据
x=iris.data
y=iris.target
x1=iris.data[:50,:4]
x2=iris.data[50:100,:4]
x3=iris.data[100:150,:4]
#print(x1)
#print(x2)
#print(x3)
print(iris)

#求每内样本的均值向量
x1_mean=x1.mean(axis=0)
x2_mean=x2.mean(axis=0)
x3_mean=x3.mean(axis=0)
print(x1_mean)
print(x2_mean)
print(x3_mean)

#求整个数据集的均值向量
x_mean=x.mean(axis=0)
print(x_mean)
print((x1_mean+x2_mean+x3_mean)/3)

#计算全局散度矩阵
x=x-x.mean(axis=0)
st=np.zeros((4,4))
print(x)
print(st)

for i in range(len(x)):
    st+=np.dot(x[i,np.newaxis].T,x[i,np.newaxis])
print(st)
#print(st/150)
print(np.dot(x.T,x))

#计算类内散度矩阵
x1=x1-x1_mean
sw1=np.zeros((4,4))
for i in range(len(x1)):
    sw1+=np.dot(x1[i,np.newaxis].T,x1[i,np.newaxis])
x2=x2-x2_mean
sw2=np.zeros((4,4))
for i in range(len(x2)):
    sw2+=np.dot(x2[i,np.newaxis].T,x2[i,np.newaxis])
x3=x3-x3_mean
sw3=np.zeros((4,4))
for i in range(len(x3)):
    sw3+=np.dot(x3[i,np.newaxis].T,x3[i,np.newaxis])
print(sw1)
print(sw2)
print(sw3)
sw=sw1+sw2+sw3
print(sw)

#计算内间散度矩阵
sb=st-sw
print(sb)

#svd求sw的逆矩阵
U,A,V=np.linalg.svd(sw)
#求A的逆矩阵
SA=np.linalg.inv(np.diag(A))
print(SA)
#计算sw的逆
VSAU=np.dot(np.dot(V.T,SA),U.T)
print(VSAU)
#计算sw-1*sb
Swb=np.dot(VSAU,sb)
print(Swb)

#特征值分解
val,vector=np.linalg.eig(Swb)
#只保留实数部分
val=np.real(val)
vector=np.real(vector)
print(val)
print(vector)

#排序
index=np.argsort(-val)
print(index)

#假设降到二维
print(index[:2])
Z=vector[:,index[:2]]
print(Z)

#降维后矩阵
Y=np.dot(x,Z)
print(Y)

#画图
for i in range(len(y)):
    if y[i]==0:
        plt.scatter(Y[i][0],Y[i][1],c='r',marker='.')
    elif y[i]==1:
        plt.scatter(Y[i][0],Y[i][1],c='b',marker='.')
    else:
        plt.scatter(Y[i][0],Y[i][1],c='g',marker='.')
plt.show()

函数

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
import random

class LinearDiscriminantAnalysis():
    def __init__(self,n_components):
        #特征数
        self.n_components=n_components
        #类别数
        self.classes=0
        #整个数据集的均值向量
        self.X_mean=0
        #特征向量
        self.eig_vectors=0
    def fit_transform(self,X,y):
        #计算样本数量以及特征数
        self.n_samples,self.n_features=X.shape
        #print(self.n_samples)
        #计算类别数，y_index存放相对于classes的类别索引
        self.classes,y_index=np.unique(y,return_inverse=True)
        #计算整体样本的均值向量
        self.X_mean=np.zeros((1,self.n_features))
        self.X_mean=np.copy(X.mean(axis=0))
       # print(self.X_mean)
        #计算全集散度矩阵
        st=np.zeros((self.n_features,self.n_features))
        for i in range(self.n_samples):
            temp=X[i]-self.X_mean
            #st+=np.dot((X[i,np.newaxis]-self.X_mean[0,np.newaxis]).T,(X[i,np.newaxis]-self.X_mean[0,np.newaxis]))
            st+=np.dot(temp[np.newaxis,:].T,temp[np.newaxis,:])
        #print(st)
        #计算每类样本的均值向量
        means=np.zeros((len(self.classes),self.n_features))
        #i代表第几个类别，从0起
        for i in range(len(self.classes)):
            #计算每类个数
            count=0
            for j in range(self.n_samples):
                if y[j]==self.classes[i]:
                    means[i]+=X[j]
                    count+=1
            means[i]=means[i]/count
        #按类别去中心化,同时求类内散度矩阵
        #print(means)
        sw=np.zeros((self.n_features,self.n_features))
        for i in range(self.n_samples):
            temp=X[i]-means[y_index[i]]
            #sw+=np.dot((X[i,np.newaxis]-means[y_index[i],np.newaxis]).T,(X[i,np.newaxis]-means[y_index[i],np.newaxis]))
            sw+=np.dot(temp[np.newaxis,:].T,temp[np.newaxis,:])
       # print(sw)
        #计算类间散度矩阵
        sb=st-sw
        #SVD求sw的逆矩阵
        U,A,V=np.linalg.svd(sw)
        #计算对角矩阵的逆矩阵
        SA=np.linalg.inv(np.diag(A))
        #计算sw的逆
        VSAU=np.dot(np.dot(V.T,SA),U.T)
        #计算SW-1Sb
        Swb=np.dot(VSAU,sb)
       # print(Swb)
        #特征值分解
        eig_vals,eig_vectors=np.linalg.eig(Swb)
        #特征值排序
        index=np.argsort(-eig_vals)
        #选取特征向量
        if self.n_components==0:
            self.n_components=len(self.classes)-1
        self.vectors=eig_vectors[:,index[:self.n_components]]
        #返回降维后的数据
        return np.dot(X,self.vectors)      

iris=datasets.load_iris()
#打乱数据
data_size=iris.data.shape[0]
#生成索引-列表
index=[i for i in range(data_size)]
#打乱
random.shuffle(index)
#通过打算的索引重新取值，相应的标签也要打乱
iris.data=iris.data[index]
iris.target=iris.target[index]
x=iris.data
y=iris.target
lda=LinearDiscriminantAnalysis(2)
Y=lda.fit_transform(X=x,y=y)
#画图
for i in range(len(y)):
    if y[i]==0:
        plt.scatter(Y[i][0],Y[i][1],c='r',marker='.')
    elif y[i]==1:
        plt.scatter(Y[i][0],Y[i][1],c='b',marker='.')
    else:
        plt.scatter(Y[i][0],Y[i][1],c='g',marker='.')
plt.show()

可能会有相反问题，这是特征向量方向的问题。

sklearn

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

iris=datasets.load_iris()
x=iris.data
y=iris.target
lda=LinearDiscriminantAnalysis(n_components=2)
Y=lda.fit(x,y).transform(x)
#画图
for i in range(len(y)):
    if y[i]==0:
        plt.scatter(Y[i][0],Y[i][1],c='r',marker='.')
    elif y[i]==1:
        plt.scatter(Y[i][0],Y[i][1],c='b',marker='.')
    else:
        plt.scatter(Y[i][0],Y[i][1],c='g',marker='.')
plt.show()