# 一、原理

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``````

``````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_)``````

``````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``````

``````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]))
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]]
``````

``````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