别再只懂PCA了！用Python手写LDA，实战鸢尾花数据集降维与分类（附完整代码）

线性判别分析实战：从数学原理到鸢尾花分类的Python实现

在数据科学领域，降维技术是处理高维数据的必备工具。提到降维，大多数人首先想到的是主成分分析(PCA)，但今天我们要探讨一个在分类任务中表现更出色的方法——线性判别分析(LDA)。与PCA不同，LDA是一种有监督的降维方法，它不仅能降低数据维度，还能最大化类别之间的区分度。

1. LDA与PCA的核心差异

**LDA（线性判别分析）和PCA（主成分分析）**都是线性降维技术，但它们的优化目标截然不同：

• PCA寻找的是数据方差最大的方向（保留最多信息）
• LDA寻找的是能最好区分不同类别的方向（最大化分类效果）

让我们用一个简单的对比表格来说明两者的关键区别：

在实际应用中，当我们的目标是分类而非单纯的数据可视化或压缩时，LDA通常是更好的选择。它能够利用已知的类别信息，找到最能区分不同类别的投影方向。

2. LDA的数学原理剖析

LDA的核心思想可以概括为：投影后类内方差最小，类间方差最大。要实现这一目标，我们需要理解几个关键概念：

2.1 类内散度矩阵(Sw)

类内散度矩阵衡量的是同一类别内样本的离散程度，计算公式为：

# 计算类内散度矩阵的Python实现 def within_class_scatter(X1, X2): # X1和X2分别代表两个类别的样本 center1 = np.mean(X1, axis=0) center2 = np.mean(X2, axis=0) cov1 = np.dot((X1 - center1).T, (X1 - center1)) cov2 = np.dot((X2 - center2).T, (X2 - center2)) Sw = cov1 + cov2 return Sw

2.2 类间散度矩阵(Sb)

类间散度矩阵衡量的是不同类别之间的分离程度：

# 计算类间散度矩阵 def between_class_scatter(center1, center2, n1, n2): mean_diff = (center1 - center2).reshape(-1, 1) Sb = np.dot(mean_diff, mean_diff.T) return Sb

2.3 求解投影方向

LDA的目标是找到一个投影方向w，使得类间散度与类内散度的比值最大化：

$$ J(w) = \frac{w^T S_b w}{w^T S_w w} $$

这个优化问题的解对应于求解广义特征值问题：

$$ S_w^{-1} S_b w = \lambda w $$

在Python中，我们可以这样实现：

def lda(X, y): # 划分不同类别样本 classes = np.unique(y) X1 = X[y == classes[0]] X2 = X[y == classes[1]] # 计算类内散度矩阵 Sw = within_class_scatter(X1, X2) # 计算类间散度矩阵 center1 = np.mean(X1, axis=0) center2 = np.mean(X2, axis=0) Sb = between_class_scatter(center1, center2, len(X1), len(X2)) # 求解广义特征值问题 eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(Sw).dot(Sb)) # 选择最大特征值对应的特征向量 w = eigenvectors[:, np.argmax(eigenvalues)] # 投影数据 X_lda = X.dot(w) return X_lda, w

注意：在实际应用中，当类内散度矩阵Sw接近奇异时，可能需要添加一个小的正则化项来保证数值稳定性。

3. 鸢尾花数据集上的LDA实战

现在，让我们将理论应用于实践，使用经典的鸢尾花数据集来演示LDA的效果。这个数据集包含三种鸢尾花（Setosa、Versicolour和Virginica）的四个特征：萼片长度、萼片宽度、花瓣长度和花瓣宽度。

3.1 数据准备与可视化

首先，我们加载数据并进行初步分析：

import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import load_iris # 加载鸢尾花数据集 iris = load_iris() X = iris.data y = iris.target feature_names = iris.feature_names target_names = iris.target_names # 可视化原始特征 plt.figure(figsize=(12, 6)) for i in range(4): plt.subplot(2, 2, i+1) for c in range(3): plt.hist(X[y == c, i], alpha=0.5, label=target_names[c]) plt.xlabel(feature_names[i]) plt.legend() plt.tight_layout() plt.show()

从特征直方图中我们可以看到，有些特征（如花瓣长度和宽度）已经能够较好地分离Setosa类别，但对Versicolour和Virginica的区分效果有限。

3.2 实现多类LDA

标准的LDA通常针对二分类问题，但我们可以将其扩展到多类情况。对于C个类别，LDA可以找到最多C-1个判别方向。

def multiclass_lda(X, y): # 总体均值 mean_total = np.mean(X, axis=0) # 计算类内散度矩阵Sw Sw = np.zeros((X.shape[1], X.shape[1])) for c in np.unique(y): Xc = X[y == c] mean_c = np.mean(Xc, axis=0) Sw += np.dot((Xc - mean_c).T, (Xc - mean_c)) # 计算类间散度矩阵Sb Sb = np.zeros((X.shape[1], X.shape[1])) for c in np.unique(y): Xc = X[y == c] mean_c = np.mean(Xc, axis=0) n_c = Xc.shape[0] Sb += n_c * np.dot((mean_c - mean_total).reshape(-1, 1), (mean_c - mean_total).reshape(1, -1)) # 求解广义特征值问题 eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(Sw).dot(Sb)) # 选择前C-1个最大特征值对应的特征向量 idx = np.argsort(eigenvalues)[::-1] W = eigenvectors[:, idx[:len(np.unique(y))-1]] # 投影数据 X_lda = X.dot(W) return X_lda, W

3.3 应用LDA并可视化结果

现在，我们将LDA应用于鸢尾花数据集：

# 应用多类LDA X_lda, W = multiclass_lda(X, y) # 可视化LDA结果 plt.figure(figsize=(8, 6)) for c in range(3): plt.scatter(X_lda[y == c, 0], X_lda[y == c, 1], alpha=0.7, label=target_names[c]) plt.xlabel('LDA Component 1') plt.ylabel('LDA Component 2') plt.title('LDA Projection of Iris Dataset') plt.legend() plt.show()

从可视化结果可以看到，经过LDA降维后，三个类别在二维空间中得到了很好的分离，特别是Setosa与其他两个类别完全分开，Versicolour和Virginica也有明显的分离趋势。

4. LDA在分类任务中的实际应用

理解了LDA的原理并实现了基础版本后，让我们看看如何在实际分类任务中使用它。

4.1 结合机器学习分类器

LDA降维后的数据可以作为其他分类器的输入：

from sklearn.model_selection import train_test_split from sklearn.linear_model import LogisticRegression from sklearn.metrics import accuracy_score # 划分训练集和测试集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) # 在训练集上拟合LDA X_train_lda, W = multiclass_lda(X_train, y_train) # 训练分类器 clf = LogisticRegression() clf.fit(X_train_lda, y_train) # 在测试集上应用LDA和分类 X_test_lda = X_test.dot(W) y_pred = clf.predict(X_test_lda) # 评估性能 accuracy = accuracy_score(y_test, y_pred) print(f"Classification accuracy: {accuracy:.2f}")

4.2 与PCA的性能对比

为了展示LDA在分类任务中的优势，我们将其与PCA进行对比：

from sklearn.decomposition import PCA # PCA降维 pca = PCA(n_components=2) X_train_pca = pca.fit_transform(X_train) X_test_pca = pca.transform(X_test) # 在PCA特征上训练分类器 clf_pca = LogisticRegression() clf_pca.fit(X_train_pca, y_train) y_pred_pca = clf_pca.predict(X_test_pca) accuracy_pca = accuracy_score(y_test, y_pred_pca) print(f"PCA + LogisticRegression accuracy: {accuracy_pca:.2f}") print(f"LDA + LogisticRegression accuracy: {accuracy:.2f}")

在大多数情况下，LDA结合简单分类器的性能会优于PCA，特别是在类别可分性较强的数据集上。这是因为LDA在降维过程中已经考虑了类别信息，而PCA只关注数据方差。

4.3 参数调优与注意事项

在实际应用中，使用LDA时需要注意以下几点：

1. 正则化：当类内散度矩阵接近奇异时，可以添加一个小的正则化项：

   Sw += np.eye(Sw.shape[0]) * 1e-4

2. 特征缩放：虽然LDA不受特征尺度影响（因为是线性变换），但在可视化时统一尺度可能更有帮助。

3. 类别平衡：LDA假设所有类别的分布形状相似，在类别不平衡时可能需要调整。

4. 维度限制：对于C个类别，LDA最多能产生C-1个判别特征。

5. 非线性扩展：对于非线性可分数据，可以考虑核LDA等扩展方法。

5. 高级话题与扩展应用

掌握了LDA的基础后，我们可以探讨一些更高级的应用场景和变体。

5.1 增量LDA处理大数据

对于大规模数据集，我们可以使用增量方式计算散度矩阵：

def incremental_lda(X_batches, y_batches): # 初始化统计量 n_features = X_batches[0].shape[1] Sw = np.zeros((n_features, n_features)) Sb = np.zeros((n_features, n_features)) mean_total = np.zeros(n_features) n_total = 0 # 第一遍计算总体均值 for X, y in zip(X_batches, y_batches): n_total += len(X) mean_total += np.sum(X, axis=0) mean_total /= n_total # 第二遍计算Sw和Sb class_stats = {} for X, y in zip(X_batches, y_batches): unique_y = np.unique(y) for c in unique_y: Xc = X[y == c] n_c = len(Xc) mean_c = np.mean(Xc, axis=0) # 更新类内散度 Sw += np.dot((Xc - mean_c).T, (Xc - mean_c)) # 更新类间散度 if c not in class_stats: class_stats[c] = {'n': 0, 'mean': np.zeros(n_features)} class_stats[c]['n'] += n_c class_stats[c]['mean'] += np.sum(Xc, axis=0) for c in class_stats: mean_c = class_stats[c]['mean'] / class_stats[c]['n'] Sb += class_stats[c]['n'] * np.dot((mean_c - mean_total).reshape(-1, 1), (mean_c - mean_total).reshape(1, -1)) # 求解投影方向 eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(Sw).dot(Sb)) idx = np.argsort(eigenvalues)[::-1] W = eigenvectors[:, idx[:len(class_stats)-1]] return W

5.2 核LDA处理非线性数据

对于非线性可分数据，我们可以使用核技巧将LDA扩展到非线性情况：

from sklearn.metrics.pairwise import rbf_kernel class KernelLDA: def __init__(self, kernel='rbf', gamma=1.0, n_components=None): self.kernel = kernel self.gamma = gamma self.n_components = n_components def fit(self, X, y): n_samples = X.shape[0] classes = np.unique(y) n_classes = len(classes) # 计算核矩阵 K = rbf_kernel(X, gamma=self.gamma) # 计算类间核矩阵 Kb = np.zeros_like(K) overall_mean = np.mean(K, axis=0) for c in classes: idx = y == c Kc = K[idx] mean_kc = np.mean(Kc, axis=0) Kb += len(Kc) * np.outer(mean_kc - overall_mean, mean_kc - overall_mean) # 计算类内核矩阵 Kw = np.zeros_like(K) for c in classes: idx = y == c Kc = K[idx] mean_kc = np.mean(Kc, axis=0) Kw += Kc.T.dot(Kc) - len(Kc) * np.outer(mean_kc, mean_kc) # 求解广义特征值问题 eigenvalues, eigenvectors = np.linalg.eig(np.linalg.pinv(Kw).dot(Kb)) idx = np.argsort(eigenvalues.real)[::-1] self.alphas = eigenvectors[:, idx[:self.n_components or (n_classes-1)]].real self.X_train = X return self def transform(self, X): K = rbf_kernel(X, self.X_train, gamma=self.gamma) return K.dot(self.alphas)

5.3 LDA在文本分类中的应用

LDA不仅适用于数值数据，在文本分类中也有广泛应用。结合TF-IDF等文本特征提取方法，LDA可以有效地降低文本数据的维度：

from sklearn.feature_extraction.text import TfidfVectorizer from sklearn.discriminant_analysis import LinearDiscriminantAnalysis # 示例文本数据 texts = ["this is a positive review", "negative experience with this product", "highly recommend this item", "would not buy again"] labels = [1, 0, 1, 0] # 1=positive, 0=negative # 提取TF-IDF特征 vectorizer = TfidfVectorizer() X_text = vectorizer.fit_transform(texts) # 应用LDA lda = LinearDiscriminantAnalysis() X_text_lda = lda.fit_transform(X_text.toarray(), labels) # 可视化文本数据的LDA投影 plt.scatter(X_text_lda[:, 0], [0]*len(X_text_lda), c=labels) plt.yticks([]) plt.title("LDA Projection of Text Data") plt.show()