LDA原理（剖析源代码，详解）
上篇⽂章我们讲解了PCA的原理，在这⾥我们先分析⼀下PCA和LDA的区别
LDA线性判别分析也是⼀种经典的降维⽅法，LDA是⼀种监督学习的降维技术，也就是说它的数据集的每个样本是有类别输出的。这点和PCA不同。PCA是不考虑样本类别输出的⽆监督降维技术。LDA的思想可以⽤⼀句话概括，就是“投影后类内⽅差最⼩，类间⽅差最⼤”。
什么意思呢？ 我们要将数据在低维度上进⾏投影，投影后希望每⼀种类别数据的投影点尽可能的接近，⽽不同类别的数据的类别中⼼之间的距离尽可能的⼤。
可能还是有点抽象，我们先看看最简单的情况。假设我们有两类数据分别为红⾊和蓝⾊，如下图所⽰，这些数据特征是⼆维的，我们希望将这些数据投影到⼀维的⼀条直线，让每⼀种类别数据的投影点尽可能的接近，⽽红⾊和蓝⾊数据中⼼之间的距离尽可能的⼤。
左边是PCA的投影效果，右边是LDA的投影效果。可以看出PCA只是考虑了整体的情况，投影到⽅差最⼤的⽅向，⽽LDA投影的时候考虑到类内⽅差减⼩，类间⽅差增⼤，这样可以更好的进⾏区分。
针对这个数据集，如果同样选择使⽤PCA，选择⽅差最⼤的⽅向作为投影⽅向，来对数据进⾏降维。那么PCA选出的最佳投影⽅向，将是图中红⾊直线所⽰的⽅向。这样做投影确实⽅差最⼤，但是是不是有其他问题。聪明的你⼀定发现了，这样做投影之后两类数据样本将混合在⼀起，将不再线性可分，甚⾄是不可分的。这对我们来说简直就是地狱，本来线性可分的样本被我们亲⼿变得不再可分。
帅⽓英俊的你也⼀定发现了，图中还有⼀条耀眼的黄⾊直线，向这条直线做投影即能使数据降维，同时还能保证两类数据仍然是线性可分的。上⾯的这个数据集如果使⽤LDA降维，出的投影⽅向就是
黄⾊直线所在的⽅向。
这其实就是LDA的思想，或者说LDA降维的⽬标：将带有标签的数据降维，投影到低维空间同时满⾜三个条件：
尽可能多地保留数据样本的信息（即选择最⼤的特征是对应的特征向量所代表的的⽅向）。
寻使样本尽可能好分的最佳投影⽅向。
投影后使得同类样本尽可能近，不同类样本尽可能远。
两者的相同点是：
1）两者均可以对数据进⾏降维。
2）两者在降维时均使⽤了矩阵特征分解的思想。
3）两者都假设数据符合⾼斯分布【正态分布】。
两者的不同点是：
1）LDA是有监督的降维⽅法，⽽PCA是⽆监督的降维⽅法
2）LDA降维最多降到类别数**k-1**的维数，⽽PCA没有这个限制。
3）LDA除了可以⽤于降维，还可以⽤于分类。（有predict⽅法，⽽PCA没有）
4）LDA选择分类性能最好的投影⽅向，⽽PCA选择样本点投影具有最⼤⽅差的⽅向。
下⾯开始进⾏分析它是如何计算出来的？
在这⾥我们还是使⽤鸢尾花的数据进⾏分析（属性较少，更直观，容易理解！）
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn import datasets
import warnings
warnings.filterwarnings('ignore')
X,y = datasets.load_iris(True)
X[:5]
array([[5.1, 3.5, 1.4, 0.2],
[4.9, 3. , 1.4, 0.2],
[4.7, 3.2, 1.3, 0.2],
[4.6, 3.1, 1.5, 0.2],
[5. , 3.6, 1.4, 0.2]])
为了验证我们计算的准确⽆误，下⾯是对原码进⾏剖析的结果，我在源代码中进⾏了修改，打印了Sw,St,Sb的准确值，⽤来⽐较如下图是源代码，在这⾥我们主要是对solver='eigen’这种情况进⾏剖析（⽐较常⽤）
class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,
TransformerMixin):
def_solve_eigen(self, X, y, shrinkage):
Sw = variance_ # within scatter
print('------Sw为',Sw)
St = _cov(X, shrinkage)# total scatter
print('******St为',St)
Sb = St - Sw # between scatter
print('++++++Sb为',Sb)
evals, evecs = linalg.eigh(Sb, Sw)
print('------evals',evals)
)[::-1][:self._max_components]
evecs = evecs[:, np.argsort(evals)[::-1]]# sort eigenvectors
# 特征值和特征向量
lda = LinearDiscriminantAnalysis(solver='eigen',n_components=2)
X_lda = lda.fit_transform(X,y)
X_lda[:5]
------Sw为 [[0.259708 0.09086667 0.164164 0.03763333]
[0.09086667 0.11308 0.05413867 0.032056 ]
[0.164164 0.05413867 0.181484 0.041812 ]
[0.03763333 0.032056 0.041812 0.041044 ]]
******St为 [[ 0.68112222 -0.04215111 1.26582 0.51282889]
[-0.04215111 0.18871289 -0.32745867 -0.12082844]
[ 1.26582 -0.32745867 3.09550267 1.286972 ]
[ 0.51282889 -0.12082844 1.286972 0.57713289]]
++++++Sb为 [[ 0.42141422 -0.13301778 1.101656 0.47519556]
[-0.13301778 0.07563289 -0.38159733 -0.15288444]
[ 1.101656 -0.38159733 2.91401867 1.24516 ]
[ 0.47519556 -0.15288444 1.24516 0.53608889]]
------evals [-2.16757273e-14 6.63220529e-15 2.85391043e-01 3.21919292e+01]
array([[6.01716893, 7.03257409],
[5.0745834 , 5.9344564 ],
[5.43939015, 6.46102462],
[4.75589325, 6.05166375],
[6.08839432, 7.24878907]])
1、总的散度矩阵
# 协⽅差
St = np.cov(X.T,bias =1)
源代码电影讲解St
array([[ 0.68112222, -0.04215111, 1.26582 , 0.51282889],
[-0.04215111, 0.18871289, -0.32745867, -0.12082844],
[ 1.26582 , -0.32745867, 3.09550267, 1.286972 ],
[ 0.51282889, -0.12082844, 1.286972 , 0.57713289]])
！！注意这⾥为什么使⽤bias=1？（重点）答案也是剖析源代码得出结果：
def_cov(X, shrinkage=None):
shrinkage ="empirical"if shrinkage is None else shrinkage
if isinstance(shrinkage,str):
if shrinkage =='auto':
sc = StandardScaler()# standardize features
X = sc.fit_transform(X)
s = ledoit_wolf(X)[0]
# rescale
s = sc.scale_[:, np.newaxis]* s * sc.scale_[np.newaxis,:]
elif shrinkage =='empirical':
s = empirical_covariance(X)
else:
raise ValueError('unknown shrinkage parameter')
从上⾯部分代码可以看出，我们对shrinkage没有设置所以会执⾏s = empirical_covariance(X)步，好奇的我就会点进去，⼀探究竟，果然！！
if assume_centered:
covariance = np.dot(X.T, X)/ X.shape[0]
else:
covariance = np.cov(X.T, bias=1)
if covariance.ndim ==0:
covariance = np.array([[covariance]])
return covariance
bias=1乖乖浮出浮出⽔⾯，这会就会发现不会有⼩的误差，准确⾄极~~
2、类内的散度矩阵
# Scatter散点图，within（内）
Sw = np.full(shape =(4,4),fill_value=0,dtype=np.float64)
for i in range(3):
Sw += np.cov(X[y == i],rowvar =False,bias =1)
Sw/=3
Sw
array([[0.259708 , 0.09086667, 0.164164 , 0.03763333],
[0.09086667, 0.11308 , 0.05413867, 0.032056 ],
[0.164164 , 0.05413867, 0.181484 , 0.041812 ],
[0.03763333, 0.032056 , 0.041812 , 0.041044 ]])
3、计算类间的散度矩阵
# Scatter between
Sb = St - Sw
Sb
array([[ 0.42141422, -0.13301778, 1.101656 , 0.47519556],
[-0.13301778, 0.07563289, -0.38159733, -0.15288444],
[ 1.101656 , -0.38159733, 2.91401867, 1.24516 ],
[ 0.47519556, -0.15288444, 1.24516 , 0.53608889]])
# scipy这个模块下的线性代数⼦模块
from scipy import linalg
4、特征值，和特征向量
eigen,ev = linalg.eigh(Sb,Sw)
display(eigen,ev)
ev = ev[:, np.argsort(eigen)[::-1]]
ev
array([-1.84103303e-14, 1.18322589e-14, 2.85391043e-01, 3.21919292e+01])
array([[ 1.54162331, -2.82590065, 0.02434685, 0.83779794],
[-2.49358543, 1.05970269, 2.18649663, 1.55005187],
[-2.86907801, 1.01439507, -0.94138258, -2.22355955],
[ 4.58628831, 0.45101349, 2.86801283, -2.83899363]])
array([[ 0.83779794, 0.02434685, -2.82590065, 1.54162331],
[ 1.55005187, 2.18649663, 1.05970269, -2.49358543],
[-2.22355955, -0.94138258, 1.01439507, -2.86907801],
[-2.83899363, 2.86801283, 0.45101349, 4.58628831]])
这⾥的ev[:, np.argsort(eigen)[::-1]]为什么要进⾏这样的转换？
看源码！就可以知道答案啦~~
class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,
TransformerMixin):
def_solve_eigen(self, X, y, shrinkage):
Sw = variance_ # within scatter
print('------Sw为',Sw)
St = _cov(X, shrinkage)# total scatter
print('******St为',St)
Sb = St - Sw # between scatter
print('++++++Sb为',Sb)
evals, evecs = linalg.eigh(Sb, Sw)
print('------evals',evals)
)[::-1][:self._max_components]
evecs = evecs[:, np.argsort(evals)[::-1]]# sort eigenvectors
X_lda[:5]
array([[6.01716893, 7.03257409],
[5.0745834 , 5.9344564 ],
[5.43939015, 6.46102462],
[4.75589325, 6.05166375],
[6.08839432, 7.24878907]])
5、删选特征向量，进⾏矩阵运算
X.dot(ev)[:,:2]
array([[ 6.01716893, 7.03257409],
[ 5.0745834 , 5.9344564 ],
[ 5.43939015, 6.46102462],
[ 4.75589325, 6.05166375],
[ 6.08839432, 7.24878907],
[ 5.65366246, 8.20566459],
.........
可以看出，原码交给我们的真的是受益匪浅，答案准确⽆误！完美！
发表评论