前言
前面一篇文章我们介绍了LDA二分类算法,这篇文章是在上一篇文章的基础上进行推广。如果推文的公式难以看懂,建议对照着西瓜书的60页看,可能我会漏一些符号的表达意义。
算法原理
在上文我们使用LDA实现了一个二分类任务。那么数据有大于2种类别,假设为C类,这时候怎么办呢?在上文我们定义的类间散度矩阵:
就不再适用,所以我们这里引入全局散度矩阵:
这里偷了一张图,可以更好的理解这个算法。
式3.35就是我们上篇博客写的“广义瑞利商”。其中的tr()为矩阵的迹,一个n×n的对角矩阵A的主对角线(从左上方至右下方的对角线)上各个元素的总和被称为矩阵A的迹(或迹数),一般记作tr(A)。 这个优化目标实际上等价于求解多个w组合成W,那么该问题就等价于求解多个上一章的优化目标,使用相同的方法,可以求得下式:
代码实现
#coding=utf-8
import numpy as np
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
# 这是sklearn中实现的LDA,待会我们会比较自己实现的LDA和它的区别
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
# k为目标
def LDA(X, y, 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
miu = np.mean(X, axis=0)
miu_classify = {}
for label in label_:
miu1 = np.mean(X_classify[label], axis=0)
miu_classify[label] = miu1
# St = np.dot((X - mju).T, X - mju)
# 计算类内散度矩阵Sw
Sw = np.zeros((len(miu), len(miu)))
for i in label_:
Sw += np.dot((X_classify[i] - miu_classify[i]).T, X_classify[i] - miu_classify[i])
#Sb = St-Sw
# 计算类内散度矩阵Sb
Sb = np.zeros((len(miu), len(miu)))
for i in label_:
Sb += len(X_classify[i]) * np.dot((miu_classify[i] - miu).reshape(
(len(miu), 1)), (miu_classify[i] - miu).reshape((1, len(miu))))
# 计算S_w^{-1}S_b的特征值和特征矩阵
eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(Sw).dot(Sb))
sorted_indices = np.argsort(eig_vals)
# 提取前k个特征向量
topk_eig_vecs = eig_vecs[:, sorted_indices[:-k - 1:-1]]
return topk_eig_vecs
def main():
iris = load_iris()
X = iris.data
y = iris.target
W = LDA(X, y, 2)
X_new = np.dot(X, W)
plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)
plt.show()
# 和sklearn的函数对比
lda = LinearDiscriminantAnalysis(n_components=2)
lda.fit(X, y)
X_new = lda.transform(X)
plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)
plt.show()
main()
可以看到使用LDA算法成功实现了多分类数据的降维。一般来说,用到LDA算法的地方不多,降维使用PCA算法更多。大概明天我会在公众号推PCA算法的原理及其实现。