KNN及降维预处理方法LDA｜PCA｜MDS

基本原理

模型介绍

knn是一种监督式学习方法

其工作机制十分简单：给定测试样本，基于某种距离度量找出训练集中与其最靠近的k个样本，然后基于这k个临近点的信息来进行该测试样本所属类别的预测。

通常而言，在分类任务中一般使用“投票法”，即选择k个样本中出现最多的类别标记作为输出结果；

在回归任务中一般使用“平均法”，即将这k个样本的实际值的平均值作为输出结果

模型分析

作为一种监督式学习模型，knn不存在显式的模型训练过程

也就是说，在训练阶段我们只需要把样本保存起来，训练时间开销几乎为0

在机器学习中，这种模式又被称为"懒惰学习"

python代码实现

python sklearn neighbors包下包含了全部knn及其拓展算法的实现

KNeighborsClassifier中比较重要的超参数有：

1. n_neighbors 表示k值
2. weights “uniform”|“distance” 表示进行投票时每个数据点的权值
3. metric “minkowski”
4. p 2 二者共同定义距离计算公式式欧氏距离，即闵氏距离p=2的情况

import numpy as np
import pandas as pd
import warnings 
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.neighbors import KNeighborsClassifierwarnings.filterwarnings("ignore")
iris_x, iris_y = load_iris(return_X_y=True)
# 为了可视化方便，我们只取iris中的前两个属性
iris_x = iris_x[:,:2]
x_train, x_test, y_train, y_test \= train_test_split(iris_x,iris_y,test_size=0.2,random_state=42)
knn_clf = KNeighborsClassifier(n_neighbors=7)
knn_clf.fit(x_train,y_train)
y_pred = knn_clf.predict(x_test)

import matplotlib as mpl
from matplotlib import pyplot as plt%matplotlib inline
mpl.rcParams["font.sans-serif"] = "SimHei"
fig,ax = plt.subplots(1,3,figsize=(12,4),facecolor="whitesmoke",edgecolor="gray")
class0 = []
class1 = []
class2 = []
for idx,cls in enumerate(iris_y):if cls == 0:class0.append(iris_x[idx,:])elif cls == 1:class1.append(iris_x[idx,:])else:class2.append(iris_x[idx,:])
ax[0].set_title("original dataset class")
ax[0].scatter(np.asarray(class0)[:,0],np.asarray(class0)[:,1],marker="*",color="r",s=3,label="class0")
ax[0].scatter(np.asarray(class1)[:,0],np.asarray(class1)[:,1],marker="*",color="orange",s=3,label="class0")
ax[0].scatter(np.asarray(class2)[:,0],np.asarray(class2)[:,1],marker="*",color="purple",s=3,label="class0")
ax[0].legend(loc="upper right")
ax[0].set_ylim([2.0,4.5])
ax[0].set_xlim([4.0,8.0])class0_pred = []
class1_pred = []
class2_pred = []
for idx,cls in enumerate(y_pred):if cls == 0:class0_pred.append(x_test[idx,:])elif cls == 1:class1_pred.append(x_test[idx,:])else:class2_pred.append(x_test[idx,:])
ax[1].set_title("original test dataset class")
ax[1].scatter(np.asarray(class0_pred)[:,0],np.asarray(class0_pred)[:,1],marker="*",color="r",s=3,label="class0")
ax[1].scatter(np.asarray(class1_pred)[:,0],np.asarray(class1_pred)[:,1],marker="*",color="orange",s=3,label="class0")
ax[1].scatter(np.asarray(class2_pred)[:,0],np.asarray(class2_pred)[:,1],marker="*",color="purple",s=3,label="class0")
ax[1].legend(loc="upper right")
ax[1].set_ylim([2.0,4.5])
ax[1].set_xlim([4.0,8.0])class0_test = []
class1_test = []
class2_test = []
for idx,cls in enumerate(y_test):if cls == 0:class0_test.append(x_test[idx,:])elif cls == 1:class1_test.append(x_test[idx,:])else:class2_test.append(x_test[idx,:])
ax[2].set_title("predict test dataset class")
ax[2].scatter(np.asarray(class0_test)[:,0],np.asarray(class0_test)[:,1],marker="*",color="r",s=3,label="class0")
ax[2].scatter(np.asarray(class1_test)[:,0],np.asarray(class1_test)[:,1],marker="*",color="orange",s=3,label="class0")
ax[2].scatter(np.asarray(class2_test)[:,0],np.asarray(class2_test)[:,1],marker="*",color="purple",s=3,label="class0")
ax[2].legend(loc="upper right")
ax[2].set_ylim([2.0,4.5])
ax[2].set_xlim([4.0,8.0])

可以看出，knn存在一定的分类误差，我们可以通过调节k的数量来使模型更精确，下面计算模型精确度：

accuracy = accuracy_score(y_pred, y_test)
$ accuracy 0.7666666666666667

降维处理

维数灾难 curse of dimensionality

经过极差对所有维度进行标准化变换后，我们可以将n个样本点视为p维单位超立方体中均匀分布的n个点

Bellman研究表明，假设样本点$X_0$处于超立方体的一个顶点位置，找到他的近邻比率$r$（近邻数量占全部数据量的比值）时超立方体各边的期望边界长度为
$$E{d}_p(r)=r^{\frac{1}{p}}$$
也就是说，对于一个$p=10$的10维空间而言，要找到$r=0.01$或者$r=0.1$的近邻，$Ed(0.01)=0.01^{\frac{1}{10}}=0.63$和$Ed(0.1)=0.1^{\frac{1}{10}}=0.8$，计算结果表明取值距离跨越了总长度的$63\%$和$80\%$

与此同时，单纯减少r的量值也不能够对超立方体各边界期望长度的减少做出多大贡献，但是同时会增大预测的方差

以上推导从一个侧面证明当维度增高时，数据点在高维度空间中的分布将不再平均，它们将更多集中在空间的边缘位置，从而导致算法的开销明显增大

下面我们将介绍两种降维的思路，其一是基于线性变化的思想，另一种则是利用低维嵌入的思想

线性变换 Linear Transformation

一般来说，最简单的维度变换即为对原始高维空间进行线性变换，给定d维空间中的样本
X = { x 1 , x 1 , ⋯ , x m } \mathbf{X}=\{\mathbf{x_1},\mathbf{x_1},\cdots,\mathbf{x_m}\} X={x1​,x1​,⋯,xm​}
变换后得到$d^{\prime }<d$维空间中的样本
$$Z=W^{T}X$$
其中$W\in {\mathbb{R}}^{d\times d^{\prime }}$是变换矩阵，$Z\in {\mathbb{R}}^{d^{\prime }\times m}$是样本在新空间中的表达

接下来，我们根据变换后对于低维子空间性质的不同要求可以构造出相应的最优化问题，并对$W$施加不同约束，

最常见的线性变换包括LDA以及PCA

LDA - 线性判别分析

传统的LDA是二分类线性判别模型，属于有监督学习的范畴，放在这里是因为我们的确使用到了线性变换的模型思路

与此同时，我们可以通过OVO或是OVR的组合模式将该二分类模型扩展至多分类模型之中

在该线性变换模型中，我们将对高维空间($n$维)中的数据点投影到$n-1$维度的超平面上，投影后的点将按照类别进行区分

相应的约束条件是：投影后类内方差最小，类间方差最大

• LDA二维示意图【实际上可能有多维】，其中“+”和“-”分别代表正例和反例

首先，我们定义数据集中的主要变量

给定数据集
$$D=\{(x_i,y_i){\}}_{i=1}^m,y_i\in (0,1)$$
令$X_i,{\mu }_i,{\Sigma }_i$分别为第$i\in (0,1)$类实例的集合、均值向量和协方差矩阵

若将数据投影到直线$\omega$上，我们可以分别得出两类样本均值和协方差矩阵的投影分别为
$${\omega }^T{\mu }_0,{\omega }^T{\mu }_1,{\omega }^T{\mathbf{\Sigma }}_0\omega ,{\omega }^T{\mathbf{\Sigma }}_1\omega$$

欲使得同类样本投影点尽可能接近，即同类样本投影点的协方差尽可能小，我们需要
$${\omega }^T{\mathbf{\Sigma }}_0\omega +{\omega }^T{\mathbf{\Sigma }}_1\omega$$
尽可能小

欲使得异类样本投影点尽可能远离，即类中心之间距离尽可能大，我们需要
$$||{\omega }^T{\mu }_0-{\omega }^T{\mu }_1|{|}_2^2$$
尽可能大

同时考虑二者，我们可以通过比率的方法构造最优化问题的目标函数
$$\begin{array}{llc}J& =& \frac{||{\omega }^T{\mu }_0-{\omega }^T{\mu }_1|{|}_2^2}{{\omega }^T{\mathbf{\Sigma }}_0\omega +{\omega }^T{\mathbf{\Sigma }}_1\omega }\\ & =& \frac{{\omega }^T({\mu }_0-{\mu }_1){({\mu }_0-{\mu }_1)}^T\omega }{{\omega }^T({\mathbf{\Sigma }}_0+{\mathbf{\Sigma }}_1)\omega }\end{array}$$

为了简化公式的书写，我们构造类内散度矩阵
$$S_w={\mathbf{\Sigma }}_0+{\mathbf{\Sigma }}_1$$
以及类间散度矩阵
$$S_b=({\mu }_0-{\mu }_1){({\mu }_0-{\mu }_1)}^T$$

则最优化问题目标函数可以重写为
$$J=\frac{{\omega }^T{S}_b\omega }{{\omega }^T{S}_w\omega }$$

由于最优化问题目标函数都是$\omega$的二次项，这表明该问题与其长度无关

不失一般性，我们令分母${\omega }^T{S}_w\omega =1$，则最优化问题可以表述如下：
$$\begin{array}{llc}& {\mathrm{argmin}}_{\omega }& -{\omega }^T{S}_b\omega \\ & s.t.& {\omega }^T{S}_w\omega =1\end{array}$$

求解该类最优化问题我们将使用拉格朗日乘子法配合KKT条件，这里不再赘述，最后得到分类直线方程为：
$$\omega ={S_{\omega }}^{-1}({\mu }_0-{\mu }_1)$$

LDA python 实现

python sklearn 包中 LinearDiscriminantAnalysis 提供LDA功能

import numpy as np
import pandas as pd 
import warnings 
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.metrics import accuracy_scorewarnings.filterwarnings("ignore")
iris_x, iris_y = load_iris(return_X_y=True)
iris_x.shape$(150,4)

可以看出，原始数据中我们拥有4个属性，也就是说降维之前我们的维度p=4

# 注意，由于受到svd影响，输出降维后的矩阵实际上是（属性值-1）和（类别值-1）二者之间的最小值
lda_model = LinearDiscriminantAnalysis(solver='svd')
lda_model.fit(iris_x,iris_y)
iris_transformed = lda_model.transform(iris_x)
iris_transformed.shape$(150,2)

使用LDA进行降维后，我们的属性来到了2个。当我们评价降维对模型带来的增益时，我们常常通过对模型本身的准确度提升来间接评价降维带来的优势

PCA - 主成分分析

主成分分析所要达到的主要目的在于能够使用一个低维的超平面更恰当的表达在高维中存在的各数据点

如何评价这样的超平面呢？我们有以下两种方式：

1. 最近重构性： 样本点到这个超平面的距离足够近
2. 最大可分性： 样本点在这个超平面上的投影能尽可能分开

有趣的是，基于以上两种判别标准，我们能够分别推导出PCA相同的最优化问题

PCA最近重构性

假定样本进行了中心化，即
$$\sum _i{x}_i=0$$
同时假定投影变换后得到的新坐标系为：
{ ω 1 , ω 1 , ⋯ , ω d } \{\omega_1,\omega_1,\cdots,\omega_d\} {ω1​,ω1​,⋯,ωd​}
其中${\omega }_i$是标准正交基向量，这意味着
$$\{\begin{array}{l}||{\omega }_i|{|}_2=1\\ {{\omega }_i}^T{\omega }_j=0,(i\ne j)\end{array}$$

PCA的目标是将维度降低到$d^{\prime }<d$，那么样本点$x_i$在低位坐标系中的投影为
$$z_i=(z_{i1},z_{i2},\cdots ,z_{i{d}^{\prime }})$$
其中$z_{ij}={{\omega }_j}^T{x}_i$是$x_i$在低维坐标系下第j维的坐标

若基于$z_i$来重构$x_i$，则会得到
$${\hat{x}}_i=\sum _{j=1}^{d^{\prime }}{z}_{ij}{\omega }_j$$

那么，我们可以构造原样本点$x_i$与基于地位样本点重构的样本点${\hat{x}}_i$之间的距离
$$\begin{array}{llc}\sum _{i=1}^m||\sum _{j=1}^{d^{\prime }}{z}_{ij}{\omega }_j-x_i|{|}_2^2& =& \sum _{i=1}^m{z_i}^T{z}_i-2\sum _{i=1}^m{z_i}^T{W}^T{x}_i+const\\ & \propto & -tr(W^T\sum _{i=1}^m{z_i}^T{x}_i)\\ & \propto & -tr(W^T\sum _{i=1}^m(W^T{x}_i{)}^T{x}_i)\\ & \propto & -tr(W^T(\sum _{i=1}^m{x}_i{x_i}^T)W)\end{array}$$

其中，$W$是一个$(d,d^{\prime })$变换矩阵，根据最近重构性、也就是重构前后样本点距离相差最小原则，以上距离应当被构造成最优化问题中的最小化结构，
考虑到${\omega }_i$是标准正交基，$\sum _{i=1}^m{x}_i{x_i}^T$是协方差矩阵，因此我们构造最优化问题如下：
$$\begin{array}{ll}\underset{W}{\mathrm{argmin}}& -tr(W^{T}X{X}^{T}W)\\ s.t.& W^{T}W=I\end{array}$$

PCA最大可分性

从最大可分性出发，我们知道，投影后的样本点可以表示为$W^T{x}_i$,
若要求所有样本点在投影后能够尽可能分开，就要求投影后样本点之间的方差最大
投影哦呼样本点之间的方差为${\sum }_i{W}^T{x}_i{x_i}^{T}W$,于是我们可以构造相应的优化问题如下：
$$\begin{array}{ll}\underset{W}{\mathrm{argmax}}& tr(W^{T}X{X}^{T}W)\\ s.t.& W^{T}W=I\end{array}$$

PCA求解及说明

我们可以发现，两种方法得到的最优化问题是对偶的。

使用拉格朗日乘子法求解最优化问题，可得：
$$X{X}^T{\omega }_i={\lambda }_i{\omega }_i$$
接下来我们对$X{X}^T$进行特征值分解，讲求的的特征值排序，取前$d^{\prime }$个特征值对应的特征向量构成$W^{\ast }$，该矩阵即为线性变换矩阵

但是我们需要注意的是，PCA模型的解释性是较低的，从另一种角度来说即为PCA降维后得到的特征是偏离原特征的实际价值的，我么只能从数学角度认为PCA处理后的特征能够更清晰地区分不同类别，但是不能够说明新生成的特征具备具体的物理意义。

PCA python实现

python sklearn decomposition.PCA 包提供PCA的相关操作

需要注意的是，PCA包中实际使用的是svd替代方阵的特征值分解，这样做的目的是减少计算开销并且尽可能小的影响精度

现在我们对sklearn.decomposition.PCA的主要参数做一个介绍：

1. n_components：这个参数可以帮我们指定希望PCA降维后的特征维度数目。最常用的做法是直接指定降维到的维度数目，此时n_components是一个大于等于1的整数。当然，我们也可以指定主成分的方差和所占的最小比例阈值，让PCA类自己去根据样本特征方差来决定降维到的维度数，此时n_components是一个（0，1]之间的数。当然，我们还可以将参数设置为"mle", 此时PCA类会用MLE算法根据特征的方差分布情况自己去选择一定数量的主成分特征来降维。我们也可以用默认值，即不输入n_components，此时n_components=min(样本数，特征数)。
2. whiten：判断是否进行白化。所谓白化，就是对降维后的数据的每个特征进行归一化，让方差都为1.对于PCA降维本身来说，一般不需要白化。如果你PCA降维后有后续的数据处理动作，可以考虑白化。默认值是False，即不进行白化。
3. svd_solver：即指定奇异值分解SVD的方法，由于特征分解是奇异值分解SVD的一个特例，一般的PCA库都是基于SVD实现的。有4个可以选择的值：{‘auto’, ‘full’, ‘arpack’, ‘randomized’}。randomized一般适用于数据量大，数据维度多同时主成分数目比例又较低的PCA降维，它使用了一些加快SVD的随机算法。 full则是传统意义上的SVD，使用了scipy库对应的实现。arpack和randomized的适用场景类似，区别是randomized使用的是scikit-learn自己的SVD实现，而arpack直接使用了scipy库的sparse SVD实现。默认是auto，即PCA类会自己去在前面讲到的三种算法里面去权衡，选择一个合适的SVD算法来降维。一般来说，使用默认值就够了。

除了这些输入参数外，有两个PCA类的成员值得关注。第一个是explained_variance_，它代表降维后的各主成分的方差值。方差值越大，则说明越是重要的主成分。第二个是explained_variance_ratio_，它代表降维后的各主成分的方差值占总方差值的比例，这个比例越大，则越是重要的主成分

import numpy as np
import pandas as pd 
import warnings
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_splitwarnings.filterwarnings("ignore")
iris_x, iris_y = load_iris(return_X_y=True)
iris_x.shape$(150,4)

pca = PCA(n_components=2)
iris_de = pca.fit_transform(iris_x)
iris_de.shape$(150,2)

多维缩放 Multiple Dimensional Scaling

多维缩放是另一种降维处理的思路，我们强调缩放前后向量的某种特性保持不变，

这种认识可以理解为我们只关心高维空间中某一个低维的分布，就好像这个低维分布是这个高维分布中的一个嵌入

因此，MDS追求原始空间中样本之间的距离在低维空间中得以保存

R代码实现

set.seed(42)
# prepare data
x1 <- runif(60,-1,1)
x2 <- runif(60,-1,1)
y <- sample(c(0,1),size=60,replace=TRUE, prob=c(0.3,0.7))
data <- data.frame(Fx1=x1, Fx2=x2, Fy=y)
# separate train and test set
sampleID <- sample(x=1:60,size=60*0.3)
x_test <- data[sampleID,]
x_train <- data[-sampleID,]
par(mfrow=c(2,2),mar=c(4,6,4,4))
plot(data[,1:2],pch=data[,3]+1,cex=0.8,xlab="x1",ylab="x2",main="total sample")
plot(x_train[,1:2],pch=x_train[,3]+1,cex=0.8,xlab="x1",ylab="x2",main="train and test sample")
points(x_test[,1:2],pch=x_test[,3]+16,col=2,cex=0.8)# use all data 
err_ratio <- vector()
for(i in 1:30){knn_fit <- knn(train=data[,1:2],test=data[,1:2],cl=as.factor(data[,3]),k=i)cross_table <- table(data[,3],knn_fit)err_ratio <- c(err_ratio,(1-sum(diag(cross_table))/sum(cross_table))*100)
}
plot(err_ratio, type="l", xlab="number of K", ylab="error rate %", main="k and error rate",ylim=c(0,80),hold=TRUE)# offsetting method 旁置法
err_ratio_offset <- vector()
for(i in 1:30){knn_fit_offset <- knn(train=x_train[,1:2],test=x_test[,1:2],cl=as.factor(x_train[,3]),k=i)cross_table_offset <- table(x_test[,3],knn_fit_offset)err_ratio_offset <- c(err_ratio_offset, (1-sum(diag(cross_table_offset))/sum(cross_table_offset))*100)
}
lines(1:30,err_ratio_offset,lty=2,col=2)# cross validation 留一法
error_ratio_cv <- vector()
for(i in 1:30){knn_fit_cv <- knn.cv(train=data[,1:2],cl=as.factor(data[,3]),k=i)cross_table_cv <- table(data[,3],knn_fit_cv)error_ratio_cv <- c(error_ratio_cv,(1-sum(diag(cross_table_cv))/sum(cross_table_cv))*100)
}
lines(1:30,error_ratio_cv,lty=3,col=4)# knn regression 回归
set.seed(24)
x3 <- runif(60,-1,1)
x4 <- runif(60,-1,1)
y2 <- runif(60,10,20)
data_reg <- data.frame(Fx3=x3, Fx4=x4, Fy=y2)
sampleID <- sample(x=1:60, size=60*0.3)
x2_test <- data_reg[sampleID,]
x2_train <- data_reg[-sampleID,]mse_vector <- vector()
for(i in 1:30){knn_fit_reg <- knn(train=x2_train[,1:2],test=x2_test[,1:2],cl=x2_train[,3],k=i,prob=FALSE)# due to the fact that res of knn is a factor, # we need to change it back to numeric vectorknn_fit_reg <- as.double(as.vector(knn_fit_reg))mse <- sum((x2_test[,3]-knn_fit_reg)^2)/length(x2_test[,3])print((x2_test[,3]-knn_fit_reg)^2)print(mse)mse_vector <- c(mse_vector,mse)
}
plot(mse_vector,type="l", xlab="k", ylab="mse",main="k-mse relationship",ylim=c(0,80))