线性判别分析（lda）

LDA是一种经典的线性学习方法，在二分类问题中由Fisher（贼屌）提出，亦称“Fisher”判别分析。
LDA思想就是给定训练样例即，设法将样例投影到一条直线上，使得同类型的投影点尽可能，一类样例的投影近可能远离，在对新样本进行分类时，将其投影到同样的这条直线上，再根据投影点的位置来确定样本的类别。

协方差：变量之间的相关程度，绝对值越大相关程度越大，若是0全都相互的独立

从公式上看，协方差是两个变量与自身期望做差再相乘，然后对乘积取期望。也就是说，当其中一个变量的取值大于自身期望，另一个变量的取值也大于自身期望时，即两个变量的变化趋势相同，此时，两个变量之间的协方差取正值。反之，即其中一个变量大于自身期望时，另外一个变量小于自身期望，那么这两个变量之间的协方差取负值。
正如上图所示，当x与y变化趋势一致时，两个变量与自身期望之差同为正或同为负，其乘积必然为正，所以其协方差为正；反之，其协方差为负。所以协方差的正负性反映了两个变量的变化趋势是否一致。再者，当x和y在某些时刻变化一致，某些时刻变化不一致时，如下图所示，在第一个点，x与y虽然变化，但是y的变化幅度远不及x变化幅度大，所以其乘积必然较小，在第二个点，x与y变化一致且变化幅度都很大，因此其乘积必然较大，在第三个点，x与y变化相反，其乘积为负值，这类点将使其协方差变小，因此，我们可以认为协方差绝对值大小反映了两个变量变化的一致程度。因此，两个变量相关系数的定义为协方差与变量标准差乘积之比。

总的来说，协方差反映了两个变量之间的相关程度。
协方差矩阵就是多变量的协方差按矩阵形式写出：

三维协方差矩阵

上文中，我们解释了协方差代表了不同维度之间的相关关系，如果说某些维度之间没有相关关系，则协方差为0，那么，以2维数据为例，我们来看一下，当不同维度之间数据没有相关关系时，即协方差矩阵为单位阵时，数据分布的整体形状。

数据协方差矩阵为单位阵时，该组数据被称为白数据，白数据在很多场合都有应用，比如在数据传输加密中，将原始数据转化成白数据，切断不同维度之间的关联关系，在访问数据时，再对数据进行解密。现在我们一起来看一下，怎么将白数据转化成真实观察数据的线性变换。协方差矩阵表示了不同方向上数据分布的离散程度，现在我们试图用一个向量和一个数字来表示协方差矩阵，向量应该指向数据分布最离散的方向，而这个数应该是该方向上数据投影的方差。首先，我们假设一个向量，则数据集D在该向量上的投影可以表示为，而数据集在该方向上的投影方差为：

二分类lda

给定数据集 D = { x i , y i } i = 1 m , y i ∈ { 0 , 1 } , 令 X i , μ i , ∑ i D =\{\boldsymbol x_i,y_i\}^m_{i=1}, y_i\in\{0, 1\},令X_i, \mu_i,\sum_i D={xi​,yi​}i=1m​,yi​∈{0,1},令Xi​,μi​,∑i​分别表示第 i , i ∈ { 0 , 1 } i, i\in\{0, 1\} i,i∈{0,1}类集合，均值向量，协方差矩阵。这里要注意 x i = { x 0 , x 1 } T , \boldsymbol x_i=\{x_0,x_1\}^T, xi​={x0​,x1​}T,，若将数据，若将数据投影到直线$\mathbf{w}$上，则两类样本的中心在直线上的投影分别为$w^T{\mu }_0和w^T{\mu }_1$；若将所有样本点都投影到直线上，则两类样本的协方差分别为${\mathbf{w}}^{\mathbf{T}}{\sum }_0\mathbf{w}和w^T{\sum }_{1}w$由于直线是一维空间，因此$w^T{\mu }_0{w}^T{\mu }_1,{\mathbf{w}}^{\mathbf{T}}{\sum }_0\mathbf{w}和w^T{\sum }_{1}w$均为实数。于是同类样本让${\mathbf{w}}^{\mathbf{T}}{\sum }_0\mathbf{w}+w^T{\sum }_{1}w$尽可能小，让$||w^T{\mu }_0-w^T{\mu }_1||2$尽可能大，即同类物体更近，异类更远。
于是由最大化目标：
$$\begin{array}{rl}J& =\frac{||w^T{\mu }_0-w^T{\mu }_1||2}{{\mathbf{w}}^{\mathbf{T}}{\sum }_0\mathbf{w}+w^T{\sum }_{1}w}\\ & \\ & =\frac{w^T({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w}{{\mathbf{w}}^{\mathbf{T}}{\sum }_0\mathbf{w}+w^T{\sum }_{1}w}\end{array}$$
定义类内散度矩阵
$$\begin{array}{rl}{S}_w& =\sum _0+\sum _1\\ & =\sum _{x\in X_0}(x-{\mu }_0)(x-{\mu }_0{)}^T+\sum _{x\in X_1}(x-{\mu }_1)(x-{\mu }_1{)}^T\end{array}$$
类间散度
$$S_b=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^T$$
于是重写
$$J=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$$

这就是欲使最大化的目标，即$S_b和S_w$的广义瑞利商
注意上式分子分母都有w的二次项，因此式子的解和w的长度无关，之和他方向有关，因此令$w^T{S}_{w}w=1$，那么问题变成了
$$\begin{gathered} mi{n}_w-w^T{S}_{b}w \\ s.t.w^T{S}_{w}w=1 \end{gathered}$$根据拉格朗日乘子法，我们构建函数
$$L(w)=-w^T{S}_{b}w+\lambda (w^T{S}_{w}w-1)$$
$$\frac{\partial (L(w))}{w}=-2S_{b}w+2\lambda S_{w}w$$令倒数为0，则上面式子等价于$${\mathbf{S}}_{\mathbf{b}}\mathbf{w}=\mathbf{\lambda }{\mathbf{S}}_{\mathbf{w}}\mathbf{w}$$
其中$\lambda$是拉格朗日乘子，注意到$S_{b}w$的方向恒为${\mu }_0-{\mu }_1$不妨令
$$S_{b}w=\lambda ({\mu }_0-{\mu }_1)$$
整理公式可以得到
$$w=S_w^{-1}({\mu }_0-{\mu }_1)$$
考虑到数值的稳定性，通常我们要对$S_w$进行奇异值分解，即$S_w=U\sum V^T$再由$S_w^{-1}=V{\sum }^{-1}{U}^T$得到$S_w^{-1}$，奇异值分解一个很重要的性质就是进行降为，同时还会保留原有的性质，是个非常重要的东西

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets.samples_generator import make_classification



def LDA(X, y):
    X1 = np.array([X[i] for i in range(len(X)) if y[i] == 0])
    X2 = np.array([X[i] for i in range(len(X)) if y[i] == 1])
    len1 = len(X1); len2 = len(X2)
    mju1 = np.mean(X1, axis=0)  # 求中心点
    mju2 = np.mean(X2, axis=0)
    cov1 = np.dot((X1 - mju1).T, (X1 - mju1))
    cov2 = np.dot((X2 - mju2).T, (X2 - mju2))
    Sw = cov1 + cov2
    w = np.dot(np.mat(Sw).I, (mju1 - mju2).reshape((len(mju1), 1)))  # 计算w
    X1_new = func(X1, w)
    X2_new = func(X2, w)
    y1_new = [1 for i in range(len1)]
    y2_new = [2 for i in range(len2)]
    return X1_new, X2_new, y1_new, y2_new

def func(x, w):
    return np.dot((x), w)


if '__main__' == __name__:
    X, y = make_classification(n_samples=500, n_features=2, n_redundant=0, n_classes=2,

                               n_informative=1, n_clusters_per_class=1, class_sep=0.5, random_state=10)

    X1_new, X2_new, y1_new, y2_new = LDA(X, y)
    # plt.scatter(X[:, 0], X[:, 1], marker='o', c=y)
    # plt.show()
    plt.plot(X1_new, y1_new, 'b*')
    plt.plot(X2_new, y1_new, 'ro')
    plt.show()

多分类lda

假设存在N个类，所有样本一共m条，定义全局散度矩阵
$$\begin{array}{rl}{S}_t& =S_b+S_w\\ & =\sum _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T\end{array}$$
$\mu$是所有样本的均值，雷内散度矩阵$S_w$定义为每个类别的散度矩阵之和，即
$$S_w=\sum _{i=1}^N{S}_{w_i}$$其中
$$S_{w_i}=\sum _{x\in X_i}(x-{\mu }_i)(x-{\mu }_i{)}^T$$
$X_i,{\mu }_i$分别是第i类样本和第i类样本的均值，
结合以上式子，得到
$$\begin{array}{rl}{S}_b& =S_t-S_w\\ & =\sum _{i=1}^N{m}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T\end{array}$$
其中，$m_i$是第i类样本的数量
$$\underset{W}{\max }\frac{tr(W_T{S}_{b}W)}{tr(W_T{S}_{w}W)}$$
tr表示矩阵的迹，$W\in R^{d\ast (N-1)}$，上式子，可以用广义特征值求解
$$S_{b}W=\lambda S_{w}W$$
W的闭式解是$S_w^{-1}{S}_b$的d个最大非零广义特征值所对应的特征向量组成的矩阵，从而达到了降为的目的。

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


def LDA_dimensionality(X, y, k):
    '''
    X为数据集，y为label，k为目标维数
    '''
    label_ = list(set(y))

    X_classify = {}

    for label in label_:
        X1 = np.array([X[i] for i in range(len(X)) if y[i] == label])
        X_classify[label] = X1

    mju = np.mean(X, axis=0)
    mju_classify = {}

    for label in label_:
        mju1 = np.mean(X_classify[label], axis=0)
        mju_classify[label] = mju1

    #St = np.dot((X - mju).T, X - mju)

    Sw = np.zeros((len(mju), len(mju)))  # 计算类内散度矩阵
    for i in label_:
        Sw += np.dot((X_classify[i] - mju_classify[i]).T,
                     X_classify[i] - mju_classify[i])

    # Sb=St-Sw

    Sb = np.zeros((len(mju), len(mju)))  # 计算类内散度矩阵
    for i in label_:
        Sb += len(X_classify[i]) * np.dot((mju_classify[i] - mju).reshape(
            (len(mju), 1)), (mju_classify[i] - mju).reshape((1, len(mju))))

    eig_vals, eig_vecs = np.linalg.eig(
        np.linalg.inv(Sw).dot(Sb))  # 计算Sw-1*Sb的特征值和特征矩阵

    sorted_indices = np.argsort(eig_vals)
    topk_eig_vecs = eig_vecs[:, sorted_indices[:-k - 1:-1]]  # 提取前k个特征向量
    return topk_eig_vecs


if '__main__' == __name__:

    iris = load_iris()
    X = iris.data
    y = iris.target
    W = LDA_dimensionality(X, y, 2)
    X_new = np.dot((X), W)
    plt.figure(1)
    plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)  # 求每个簇的中心点，新样本点根据距离的远近确定类别

    # 与sklearn中的LDA函数对比
    lda = LinearDiscriminantAnalysis(n_components=2)
    lda.fit(X, y)
    X_new = lda.transform(X)
    print(X_new)
    plt.figure(2)
    plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)

    plt.show()