线性判别分析(Linear Discriminant Analysis, LDA)

主要思想

LDA将$D$维特征$\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N]\in {\mathbb{R}}^{D\times N}$（${\mathbf{x}}_i\in {\mathbb{R}}^D$）映射到$d(d\ll D)$维空间中（$\mathbf{Z}=[{\mathbf{z}}_1,{\mathbf{z}}_2,\cdots ,{\mathbf{z}}_N]={\mathbf{W}}^T\mathbf{X}\in {\mathbb{R}}^{d\times N}$, $\mathbf{W}=[{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d]\in {\mathbb{R}}^{D\times d}$，${\mathbf{w}}_i\in {\mathbb{R}}^D$），使得在降维后的空间中，尽量使同一类的数据聚集，不同类的数据尽可能分散，因此对于每个特征${\mathbf{x}}_i$还需要知道其对应的标签 y i ∈ { 1 , 2 , ⋯   , C } y_i\in\{1,2,\cdots,\mathcal{C}\} yi​∈{1,2,⋯,C}。

推导方法

此处我们不考虑二分类的情况，直接考虑多分类的情况。根据降维思想，降维后的$\mathbf{Z}$满足同类更紧凑，而异类更分散，那么衡量“紧凑”和“分散”可以定义类内散度矩阵和类间散度矩阵：
$$\begin{array}{rl}{\mathbf{S}}_w& =\sum _{c=1}^{\mathcal{C}}\frac{N_c}{N}{\mathbf{S}}_w^{(c)}\\ & =\sum _{c=1}^{\mathcal{C}}\frac{N_c}{N}\frac{1}{N_c}\sum _{i=1}^{N}1[y_i=c]({\mathbf{x}}_i-{\mu }_c)({\mathbf{x}}_i-{\mu }_c{)}^T\\ {\mathbf{S}}_b& =\frac{1}{2}\sum _{c_1,c_2=1}^{\mathcal{C}}\frac{N_{c_1}}{N}\frac{N_{c_2}}{N}{\mathbf{S}}_b^{(c_1,c_2)}\\ & =\frac{1}{2}\sum _{c_1,c_2=1}^{\mathcal{C}}\frac{N_{c_1}}{N}\frac{N_{c_2}}{N}({\mu }_{c_1}-{\mu }_{c_2})({\mu }_{c_1}-{\mu }_{c_2}{)}^T\\ & =\frac{1}{2}\sum _{c_1=1}^{\mathcal{C}}\frac{N_{c_1}}{N}\sum _{c_2=1}^{\mathcal{C}}\frac{N_{c_2}}{N}({\mu }_{c_1}-{\mu }_{c_2})({\mu }_{c_1}-{\mu }_{c_2}{)}^T\\ & =\frac{1}{2}\left(\sum _{c_1=1}^{\mathcal{C}}\frac{N_{c_1}}{N}{\mu }_{c_1}{\mu }_{c_1}^T-\mu {\mu }^T-\mu {\mu }^T+\sum _{c_2=1}^{\mathcal{C}}\frac{N_{c_2}}{N}{\mu }_{c_2}{\mu }_{c_2}^T\right)\\ & =\sum _{c_1=1}^{\mathcal{C}}\frac{N_{c_1}}{N}{\mu }_{c_1}{\mu }_{c_1}^T-\mu {\mu }^T\\ & =\sum _{c_1=1}^{\mathcal{C}}\frac{N_{c_1}}{N}({\mu }_{c_1}{\mu }_{c_1}^T-\mu {\mu }^T+2\mu {\mu }^T-{\mu }_{c_1}{\mu }^T-\mu {\mu }_{c_1}^T)\\ & =\sum _{c_1=1}^{\mathcal{C}}\frac{N_{c_1}}{N}({\mu }_{c_1}{\mu }_{c_1}^T+\mu {\mu }^T-{\mu }_{c_1}{\mu }^T-\mu {\mu }_{c_1}^T)\\ & =\sum _{c=1}^{\mathcal{C}}\frac{N_c}{N}\left({\mu }_c-\mu \right){\left({\mu }_c-\mu \right)}^T\end{array}$$
其中
$$\begin{array}{rl}{N}_c& =\sum _{i=1}^{N}1[y_i=c]\\ {\mu }_c& =\frac{1}{N_c}\sum _{i=1}^{N}1[y_i=c]{\mathbf{x}}_i\\ \mu & =\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_i\\ & =\sum _{c=1}^{\mathcal{C}}\frac{N_c}{N}{\mu }_c\end{array}$$
以上的${\mathbf{S}}_b$和${\mathbf{S}}_w$计算的是$\mathbf{X}$，而目标是找到$\mathbf{W}$使得$\mathbf{Z}={\mathbf{W}}^T\mathbf{X}$在同类中更加紧凑，在异类中更分散。优化目标有很多种表示方法：

求和优化目标

$$\begin{gathered} \underset{\mathbf{W}}{\mathrm{argmax}}\frac{trace\left({\mathbf{W}}^T{\mathbf{S}}_b\mathbf{W}\right)}{trace\left({\mathbf{W}}^T{\mathbf{S}}_w\mathbf{W}\right)} \\ \iff \{\begin{array}{c}\underset{\mathbf{W}}{\mathrm{argmax}}trace\left({\mathbf{W}}^T{\mathbf{S}}_b\mathbf{W}\right)\\ s.t.{\mathbf{W}}^T{\mathbf{S}}_w\mathbf{W}=\mathbb{I}\end{array} \\ \iff \{\begin{array}{c}\underset{\mathbf{W}}{\mathrm{argmax}}trace\left({\mathbf{V}}^T{\mathbf{S}}_w^{-1/2}{\mathbf{S}}_b{\mathbf{S}}_w^{-1/2}\mathbf{V}\right)\\ s.t.{\mathbf{V}}^T\mathbf{V}=\mathbb{I}\\ \mathbf{V}={\mathbf{S}}_w^{1/2}\mathbf{W}\end{array} \end{gathered}$$
那么以上的求解与PCA所用的方法一致了，然后再$\mathbf{W}={\mathbf{S}}_w^{-1/2}\mathbf{V}$就OK啦！需要说明的是，由于$rank({\mathbf{S}}_b)\le \mathcal{C}-1$，所以最终求得的$\mathbf{V}$最多就$\mathcal{C}-1$个特征向量。

求积优化目标

$$\begin{gathered} \underset{\mathbf{W}}{\mathrm{argmax}}\prod _{i=1}^d\frac{{\mathbf{w}}_i^T{\mathbf{S}}_b{\mathbf{w}}_i}{{\mathbf{w}}_i^T{\mathbf{S}}_w{\mathbf{w}}_i} \\ \iff \{\begin{array}{c}\underset{\mathbf{W}}{\mathrm{argmax}}{\prod }_{i=1}^d{\mathbf{v}}_i^T{\mathbf{S}}_w^{-1/2}{\mathbf{S}}_b{\mathbf{S}}_w^{-1/2}{\mathbf{v}}_i\\ s.t.{\mathbf{V}}^T\mathbf{V}=\mathbb{I}\\ \mathbf{V}={\mathbf{S}}_w^{1/2}\mathbf{W}\end{array} \end{gathered}$$
更简单的一种是${\mathbf{S}}_w^{-1/2}{\mathbf{S}}_b{\mathbf{S}}_w^{-1/2}$求出的特征向量$\mathbf{v}$不用左乘${\mathbf{S}}_w^{-1/2}$得到$\mathbf{w}$，可以直接通过求解${\mathbf{S}}_w^{-1}{\mathbf{S}}_b$得到$\mathbf{w}$。

from stat import S_IFBLK
import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

def plot_2d(X_r, y, target_names, name):
    if 'PCA' in name:
        plt.subplot(1,3,1)
    else:
        if 'Sklearn' in name:
            plt.subplot(1,3,2)
        else:
            plt.subplot(1,3,3)
    
    colors = ["navy", "turquoise", "darkorange"]
    lw = 2

    for color, i, target_name in zip(colors, [0, 1, 2], target_names):
        plt.scatter(
            X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=0.8, lw=lw, label=target_name
        )
    plt.legend(loc="best", shadow=False, scatterpoints=1)
    plt.xlabel(f"{name} of IRIS dataset")

def sklearn_lda(X, y, target_names):
    lda = LinearDiscriminantAnalysis(n_components=2)
    X_r2 = lda.fit(X, y).transform(X)
    plot_2d(X_r2, y, target_names, 'Sklearn LDA')
    return X_r2

def my_lda(X, y, target_names):
    Sw = 0.0
    Sb = 0.0
    mu = np.mean(X, axis=0)
    for i in range(3):
        X_c  = X[y==i,:]
        mu_c = np.mean(X_c, axis=0)
        N_c  = X_c.shape[0]

        Sw = Sw + (X_c-mu_c).T @ (X_c-mu_c)

        Sb = Sb + N_c * np.reshape(mu_c - mu, [-1, 1]) @ np.reshape(mu_c - mu, [1, -1])
    S = np.linalg.inv(Sw) @ Sb

    values, vectors = np.linalg.eig(S)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,3)]]
    
    X_r = np.matmul(X-mu, W)
    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致
    X_r[:,0] = - X_r[:,0] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. LDA')

    return X_r

def my_pca(X, y, target_names):
    n_components=2
    # 去中心化
    N = X.shape[0]
    X_mean = np.mean(X, axis=0)
    X_ = X - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. PCA')

    return X_r



if __name__ == '__main__':
    iris = datasets.load_iris()

    X = iris.data
    y = iris.target
    target_names = iris.target_names

    plt.figure()

    

    X_r2 = my_pca(X, y, target_names)
    X_r1 = my_lda(X, y, target_names)
    X_r1 = sklearn_lda(X, y, target_names)

    plt.show()

核技巧

之前的问题实质上是求解
$${\mathbf{S}}_b{\mathbf{w}}_k={\lambda }_k{\mathbf{S}}_w{\mathbf{w}}_k$$
的广义特征向量。首先将核化后的${\mathbf{S}}_b$和${\mathbf{S}}_w$写出
$$\begin{array}{rl}{\mathbf{S}}_w^{\varphi }& =\frac{1}{N}\sum _{c=1}^{\mathcal{C}}\sum _{i=1}^{N}1\left[y_i=c\right]\left(\varphi \left({\mathbf{x}}_i\right)-\frac{1}{N_c}{\mathbf{X}}_{\varphi }{1}_{N\times 1}^{\left(c\right)}\right){\left(\varphi \left({\mathbf{x}}_i\right)-\frac{1}{N_c}{\mathbf{X}}_{\varphi }{1}_{N\times 1}^{\left(c\right)}\right)}^T\\ & =\frac{1}{N}\sum _{c=1}^{\mathcal{C}}\sum _{i=1}^{N}1\left[y_i=c\right]\left(\varphi \left({\mathbf{x}}_i\right)\varphi {\left({\mathbf{x}}_i\right)}^T-{\mathbf{X}}_{\varphi }\frac{1_{N\times 1}^{\left(c\right)}}{N_c}\varphi {\left({\mathbf{x}}_i\right)}^T-\varphi \left({\mathbf{x}}_i\right)\frac{1_{1\times N}^{\left(c\right)}}{N_c}{\mathbf{X}}_{\varphi }^T+{\mathbf{X}}_{\varphi }\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c^2}{\mathbf{X}}_{\varphi }^T\right)\\ & =\frac{1}{N}\left({\mathbf{X}}_{\varphi }{\mathbf{X}}_{\varphi }^T-{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T-{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T+{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T\right)\\ & =\frac{1}{N}{\mathbf{X}}_{\varphi }\left(\mathbb{I}-\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T\end{array}$$
$$\begin{array}{rl}{\mathbf{S}}_b^{\varphi }& =\frac{1}{N}\sum _{c=1}^{\mathcal{C}}{N}_c\left(\frac{1}{N_c}{\mathbf{X}}_{\varphi }{1}_{N\times 1}^{\left(c\right)}-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}\right){\left(\frac{1}{N_c}{\mathbf{X}}_{\varphi }{1}_{N\times 1}^{\left(c\right)}-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}\right)}^T\\ & =\frac{1}{N}\left({\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T-{\mathbf{X}}_{\varphi }\left(\frac{1_{N\times 1}}{N}\sum _{c=1}^{\mathcal{C}}{1}_{1\times N}^{\left(c\right)}\right){\mathbf{X}}_{\varphi }^T-{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}{1}_{N\times 1}^{\left(c\right)}\frac{1_{1\times N}}{N}\right){\mathbf{X}}_{\varphi }^T+{\mathbf{X}}_{\varphi }\left(1_{N\times 1}{1}_{1\times N}\sum _{c=1}^{\mathcal{C}}\frac{N_c}{N^2}\right){\mathbf{X}}_{\varphi }^T\right)\\ & =\frac{1}{N}{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times 1}^{\left(c\right)}{1}_{1\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times 1}{1}_{1\times N}}{N}-\frac{1_{N\times 1}{1}_{1\times N}}{N}+\frac{1_{N\times 1}{1}_{1\times N}}{N}\right){\mathbf{X}}_{\varphi }^T\\ & =\frac{1}{N}{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right){\mathbf{X}}_{\varphi }^T\end{array}$$
而根据核技巧有
$$\begin{array}{rl}{\mathbf{w}}_k& ={\mathbf{X}}_{\varphi }{\alpha }_k\\ & =\left[\varphi \left({\mathbf{x}}_1\right),\varphi \left({\mathbf{x}}_2\right),\cdots ,\varphi \left({\mathbf{x}}_N\right)\right]\left[\begin{array}{c}{\alpha }_{1k}\\ {\alpha }_{2k}\\ ⋮\\ {\alpha }_{Nk}\end{array}\right]\end{array}$$
那么
$$\begin{gathered} \frac{1}{N}{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right){\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k={\lambda }_k\frac{1}{N}{\mathbf{X}}_{\varphi }\left(\mathbb{I}-\sum _{c=1}^{\mathcal{C}}\frac{1{}_{N\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k \\ {\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right){\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k={\lambda }_k{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }\left(\mathbb{I}-\sum _{c=1}^{\mathcal{C}}\frac{1{}_{N\times N}^{\left(c\right)}}{N_c}\right){\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k \\ \mathbf{K}\left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right)\mathbf{K}{\alpha }_k={\lambda }_k\mathbf{K}\left(\mathbb{I}-\sum _{c=1}^{\mathcal{C}}\frac{1{}_{N\times N}^{\left(c\right)}}{N_c}\right)\mathbf{K}{\alpha }_k \\ \left(\sum _{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right)\mathbf{K}{\alpha }_k={\lambda }_k\left(\mathbb{I}-\sum _{c=1}^{\mathcal{C}}\frac{1{}_{N\times N}^{\left(c\right)}}{N_c}\right)\mathbf{K}{\alpha }_k \end{gathered}$$
所以${\alpha }_k$就是${\left[\left(\mathbb{I}-{\sum }_{c=1}^{\mathcal{C}}\frac{1{}_{N\times N}^{\left(c\right)}}{N_c}\right)\mathbf{K}\right]}^{-1}\left({\sum }_{c=1}^{\mathcal{C}}\frac{1_{N\times N}^{\left(c\right)}}{N_c}-\frac{1_{N\times N}}{N}\right)\mathbf{K}$的特征向量。对一个新的数据${\mathbf{x}}_{new}$，降维到${\mathbf{w}}_k$的维度
$$\begin{array}{rl}{\mathbf{w}}_k^T\left(\varphi \left({\mathbf{x}}_{new}\right)-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}\right)& ={\alpha }_k^T{\mathbf{X}}_{\varphi }^T\varphi \left({\mathbf{x}}_{new}\right)-\frac{1}{N}{\alpha }_k^T{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{1}_{N\times 1}\\ & ={\alpha }_k^T\mathbf{K}\left(\cdot ,{\mathbf{x}}_{new}\right)-\frac{1}{N}{\alpha }_k^T\mathbf{K}{1}_{N\times 1}\end{array}$$

from stat import S_IFBLK
import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

def plot_2d(X_r, y, target_names, name):
    if 'PCA' in name:
        plt.subplot(1,3,1)
    else:
        if 'KLDA' in name:
            plt.subplot(1,3,3)
        else:
            plt.subplot(1,3,2)
    
    colors = ["navy", "turquoise", "darkorange"]
    lw = 2

    for color, i, target_name in zip(colors, [0, 1, 2], target_names):
        plt.scatter(
            X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=0.8, lw=lw, label=target_name
        )
    plt.legend(loc="best", shadow=False, scatterpoints=1)
    plt.xlabel(f"{name} of IRIS dataset")

def my_klda(X, y, target_names):
    gamma = 10.0
    # exp(-||x_i-x_j||^2*gamma)
    N = X.shape[0]
    K0 = np.sum(X**2, axis=1, keepdims=True)
    K0 = np.exp(-(K0 + K0.T - 2 * np.matmul(X, X.T))*gamma)

    Ic = 0
    for i in range(3):
        X_c  = X[y==i,:]
        N_c  = X_c.shape[0]
        lc = np.zeros([N, 1])
        lc[y==i,:] = 1.0
        Ic = Ic + lc @ lc.T / N_c
    
    Sb_ = (Ic - np.ones([N, N]) / N) @ K0
    Sw_ = (np.eye(N) - Ic) @ K0

    S = np.linalg.inv(Sw_) @ Sb_

    values, vectors = np.linalg.eig(S)
    idxs = np.argsort(values)
    
    alphas = vectors[:,[idxs[-i] for i in range(1,3)]]
    
    X_r = alphas.T @ K0 - 1. / N * alphas.T @ K0 @ np.ones([N, 1])
    X_r = X_r.T

    plot_2d(X_r, y, target_names, 'My Imple. KLDA')

    return X_r
    

    


def my_lda(X, y, target_names):
    Sw = 0.0
    Sb = 0.0
    mu = np.mean(X, axis=0)
    for i in range(3):
        X_c  = X[y==i,:]
        mu_c = np.mean(X_c, axis=0)
        N_c  = X_c.shape[0]

        Sw = Sw + (X_c-mu_c).T @ (X_c-mu_c)

        Sb = Sb + N_c * np.reshape(mu_c - mu, [-1, 1]) @ np.reshape(mu_c - mu, [1, -1])
    S = np.linalg.inv(Sw) @ Sb

    values, vectors = np.linalg.eig(S)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,3)]]
    
    X_r = np.matmul(X-mu, W)
    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致
    X_r[:,0] = - X_r[:,0] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. LDA')

    return X_r

def my_pca(X, y, target_names):
    n_components=2
    # 去中心化
    N = X.shape[0]
    X_mean = np.mean(X, axis=0)
    X_ = X - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. PCA')

    return X_r



if __name__ == '__main__':
    iris = datasets.load_iris()

    X = iris.data
    y = iris.target
    target_names = iris.target_names

    plt.figure()

    

    X_r2 = my_pca(X, y, target_names)
    X_r1 = my_lda(X, y, target_names)
    X_r1 = my_klda(X, y, target_names)

    plt.show()