线性判别分析LDA用于降维
一、原理
LDA线性判别分析，LDA是一种监督学习的降维技术，LDA的思想可以用一句话概括，就是“投影后类内方差最小，类间方差最大”。
二、实现
1、基于奇异值分解svd的LDA降维
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
if __name__ == '__main__':
X_ = np.array([[-1, -1],
[-2, -1],
[-3, -2],
[1, 1],
[2, 1],
[3, 2]])
y_ = np.array([1, 1, 1, 2, 2, 2])
clf = LinearDiscriminantAnalysis(solver='svd')
clf.fit(X_, y_)
print(clf.transform(X_))
# [[-1.732], [-1.732], [-3.464], [1.732], [1.732], [3.464]] --- svd分解后代码（python）
import numpy as np
from scipy import linalg
if __name__ == '__main__':
print('Hello world!')
X_ = np.array([[-1, -1],
[-2, -1],
[-3, -2],
[1, 1],
[2, 1],
[3, 2]])
y_ = np.array([1, 1, 1, 2, 2, 2])
# class mean
classes, y = np.unique(y_, return_inverse=True) # unique -- classes list, indexes
cnt = np.bincount(y) # cnt -- cnt for every class
means = np.zeros(shape=(len(classes), X_.shape[1])) # n_c * n_f
np.add.at(means, y, X_) # every class data sum
means /= cnt[:, None] # every class means
n_c = len(classes) # class number
n_s = X_.shape[0] # samples number
priors = np.bincount(y) / float(len(y_)) # 1*n_c
xbar = np.dot(priors, means) # 1*nc * nc*nf = 1*nf class-weight means
# every class center simple
Xc = []
for idx, group in enumerate(classes):
Xg = X_[y_ == group, :]
Xc.append(Xg - means[idx])
Xc = np.concatenate(Xc, axis=0) # ns * nf center simple
std = Xc.std(axis=0) # 1*nf every feature std
std[std == 0] = 1. # ignore zero
fac = 1. / (n_s - n_c) # fac
X = np.sqrt(fac) * (Xc / std) # X standardization ns * nf
U, S, Vt = linalg.svd(X, full_matrices=False) # svd ---- [1.67303261 0.44828774] U*S*Vt=X
rank = np.sum(S > 1e-4) # feature value > 0 cnt
scaling = (Vt[:rank] / std).T / S[:rank] # Vt/std/S ---- scale
fac2 = np.sqrt((n_s * priors) * fac) # [0.8660254 0.8660254] 1*nc
X = np.dot((fac2*(means-xbar).T).T, scaling) #
_, S, Vt = linalg.svd(X, full_matrices=False) # svd ---- [2.8284e+00 8.3677e-18]
rank = np.sum(S > 1e-4 * S[0]) # feature value > 0 cnt
scaling_ = np.dot(scaling, Vt.T[:, :rank]) # Scale * Vt
X_new = np.dot(X_-xbar, scaling_)分解后代码（matlab）
clear all
clc
% 基于奇异值分解的LDA降维实现
X = [[-1, -1]; [-2, -1]; [-3, -2]; [1, 1]; [2, 1]; [3, 2]];
y = [1, 1, 1, 2, 2, 2];
labels = unique(y); % 所有类别
nc = length(labels); % 类别数量
[ns, nf] = size(X); % 样本数 特征维度
mu = zeros(nc, nf); % 每个类别的均值 nc*nf
p = zeros(1, nc); % 各类别样本所占比例
Xc = zeros(ns, nf); % 各类去中心化的数据
ind = 1; % 下标索引
for i = 1:nc
index = find(y==labels(i)); % 类别i的样本下标
ni = length(index); % 类别i的样本数
Xi = X(index, :); % 类别i的样本 n*n_f
mui = mean(Xi); % 类别i的样本均值
mu(i, :) = mui; % 类别i的均值
mui = repmat(mui, ni, 1); % 类别i的均值矩阵 ni*n_f
Xi = Xi - mui; % 去中心化 ni*n_f
Xc(ind:ind+ni-1, :) = Xi; % 各类去中心化的数据
ind = ind + ni; % 更新下标
p(i) = ni / ns; % 各类别样本所占比例
end
xbar = p * mu; % class-weight特征均值 1*nf
s = std(Xc, 1); % s=((1/n)*sum(xi-mu)^2)^(1/2) 1*nf
s(s==0) = 1; % ignore zero
fac = 1 / (ns - nc); % factor
ss = repmat(s, ns, 1); % ns*nf
Xs = sqrt(fac) * (Xc./ss); % Xs ns*nf
[U, S, V] = svd(Xs, 0); % SVD U*S*V'=Xs
V = V'; % 转置
S = diag(S)'; % 对角元素
d = sum(S>0.0001); % 特征值>0的个数
V = V(1:d, :); % 前d行
ss = repmat(s, d, 1); % d*nf
sss = repmat(S(1:d), nf, 1); % nf*d
scale = (V ./ ss)' ./ sss; % nc*nf
fac2 = sqrt(fac*ns*p); % 1*nc
fac2s = repmat(fac2, nf, 1); % nf*nc
xbars = repmat(xbar, nc, 1); % nc*nf
fac3 = (fac2s .* (mu - xbars)')'; % nc*nf
Xt = fac3 * scale; % nc*nf
[U1, S1, V1] = svd(Xt, 0);
S1 = diag(S1)'; % 对角元素
d1 = sum(S1>0.0001*S1(1)); % 特征值>0的个数
scale2 = scale * V1(:, 1:d1); % nc*d1
xbars2 = repmat(xbar, ns, 1);
Xnew = (X-xbars2) * scale2;
2、基于特征分解eig的LDA降维实现
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
if __name__ == '__main__':
X_ = np.array([[-1, -1],
[-2, -1],
[-3, -2],
[1, 1],
[2, 1],
[3, 2]])
y_ = np.array([1, 1, 1, 2, 2, 2])
clf = LinearDiscriminantAnalysis(solver='eigen')
clf.fit(X_, y_)
print(clf.transform(X_))
# [[2.121], [2.121], [4.243], [-2.121], [-2.121], [-4.243]] --- eigen分解后代码（python）
import numpy as np
from scipy import linalg
def class_means_(X, y):
classes, y = np.unique(y, return_inverse=True)
cnt = np.bincount(y)
means = np.zeros(shape=(len(classes), X.shape[1]))
np.add.at(means, y, X)
means /= cnt[:, None]
return means
def cov_(X):
X = np.asarray(X)
covariance = np.cov(X.T, bias=True) # by N-1 if bias=False; by N if bias=True
return covariance
def class_cov_(X, y, priors):
cls = np.unique(y)
cov = np.zeros(shape=(X.shape[1], X.shape[1]))
for idx, group in enumerate(cls):
Xg = X[y == group, :]
cov += priors[idx] * np.atleast_2d(cov_(Xg))
return cov
if __name__ == '__main__':
print('Hello world!')
X_ = np.array([[-1, -1],[-2, -1],[-3, -2],[1, 1],[2, 1],[3, 2]])
y_ = np.array([1, 1, 1, 2, 2, 2])
classes, y_t = np.unique(y_, return_inverse=True)
n_c = len(classes) # class number
n_s = X_.shape[0] # samples number
n_f = X_.shape[1] # feature number
max_components = min(n_c - 1, n_f) # nax dim
priors = np.bincount(y_t) / float(len(y_)) # 1*n_c
means_ = class_means_(X_, y_) # class means n_c*n_f
covariance_ = class_cov_(X_, y_, priors) # class cov n_f*nf
Sw = covariance_ # within scatter [[0.6667, 0.3333], [0.3333, 0.2222]]
St = cov_(X_) # total scatter
Sb = St - Sw # between scatter [[4.0000, 2.6667], [2.6667, 1.7778]]
evals, evecs = linalg.eigh(Sb, Sw) # only used by Sym Matrix Sb*v=lambda*Sw*v v.T*Sb*v=lambda v.T*Sw*v=1
# [0, 8] [[2.4495, 0], [-3.6742, -2.1213]]
evals2, evecs2 = linalg.eig(Sb, Sw) # any matrix [0+0.j, 8+0.j] [[-0.5547, 0], [-0.8321, 1]]
evecs = evecs[:, np.argsort(evals)[::-1]] # sort eigenvectors
X_new = np.dot(X_, evecs) # X*V n_s*n_f * n_f*n_f = n_s*n_f
print(X_new[:, :max_components]) # decrease dim
# [[2.121], [2.121], [4.243], [-2.121], [-2.121], [-4.243]]
分解后代码（matlab）
clear all
clc
% 基于特征分解的LDA降维实现
X = [[-1, -1]; [-2, -1]; [-3, -2]; [1, 1]; [2, 1]; [3, 2]];
y = [1, 1, 1, 2, 2, 2];
labels = unique(y); % 所有类别
nc = length(labels); % 类别数量
nf = size(X, 2); % 特征维度
ns = size(X, 1); % 样本数
m = mean(X); % 所有样本均值 1*nf
ms = repmat(m, ns, 1); % ns*nf
Xc = X - ms; % 去中心
St = (1/ns)*(Xc'*Xc); % 整体散度
Sw = zeros(nf, nf); % 类内散度矩阵
Sb = zeros(nf, nf); % 类间散度
for i = 1:nc
ind = find(y==labels(i)); % 类别i的样本下标
ni = length(ind); % 类别i的样本数ni
Xi = X(ind, :); % 类别i的样本 ni*nf
mi = mean(Xi); % 类别i的样本均值 1*nf
mis = repmat(mi, ni, 1); % 类别i的均值矩阵 ni*nf
Xi = Xi - mis; % 去中心化 ni*nf
Sigma = (1/ni)*(Xi'*Xi); % 协方差矩阵 nf*nf
Sw = Sw + (ni/ns)*Sigma; % 类内散度，越小越好
Beta = (ni/ns)*(mi-m)'*(mi-m); % 协方差矩阵 nf*nf
Sb = Sb + Beta; % 类间散度，越大越好
end
Sb2 = St - Sw; % 类间散度--等价Sb
[Vec, Val] = eig(Sb, Sw); % 广义特征分解
index = nf:-1:1; % 逆序
Vec = Vec(:, index); % 逆序
d = min(nc-1, nf); % 降维后的维度
R = X*Vec(:,1:d); % 降维后数据
本图文内容来源于网友网络收集整理提供，作为学习参考使用，版权属于原作者。
THE END
二维码