第11章 降维 笔记
去掉数据集中关联性不大和冗余的数据，确保不出现过度适应的前提下降低计算的成本，需要对特征进行无损规约，数学上叫降维。广泛用于模式识别、文本检索以及机器学习领域，主要分为两类，特征提取和特征筛选，前者是高维数据投影到低维空间，后者是特征子集代替原始特征集，包括特征分级和特征筛选，分级是找到优化后的特征子集。特征提取可以分成线性抽取和非线性抽取两种方法，前者是试图找到一个仿射空间能够最好的说明数据分布的变化，后者对高维非线性曲线平面分布的数据非常有效。线性特征的抽取方法：
- PCA:主成分分析，将数据映射到低维空间，使低维数据最大限度保持方差信息。
- MDS:多维尺度分析，帮助用户掌握对象之间的相对距离（模式距离），并获得低维空间的表达形式，如果以数据的协方差作为距离，可以将PCA看作MDS的一种简单的形式。
- SVD: 奇异值分解，去掉线性代数角度存在线性相关冗余的特征，PCA也是它的特例。
- ISOMAP：可以看作MDS的扩展，使用几何距离作距离，是指两点之间的最短路径。
- LLE: 局部线性嵌入法，采用局部PCA和全局特征分解，相比ISOMAP，LLE主要处理局部解，目标是为每个特征选择合适的特征子集，而前者面向处理全部的特征。
11.2 FSelector完成特征筛选
首先设定一些标准，然后挑选出满足标准的特征。
install.packages("FSelector")
library(FSelector)
library(modeldata)
data(mlc_churn)
#data(iris)
churnTrain <- mlc_churn
ind <- sample(2,nrow(churnTrain),replace = TRUE,
prob = c(0.7,0.3))
trainset <- churnTrain[ind==1,]
testset <- churnTrain[ind==2,]
#churn.bagging <- bagging(Species~., data = trainset, coob=TRUE)
weights <- random.forest.importance(churn~., trainset, importance.type = 1)
print(weights)
attr_importance
state 1.6873579
account_length -2.0627558
area_code 0.5575704
international_plan 80.7366154
...
subset<- cutoff.k(weights,5)
f <- as.simple.formula(subset,"Class");f
Class ~ number_customer_service_calls + international_plan +
total_day_minutes + total_day_charge + total_intl_calls
<environment: 0x7f86c888bcc8>
evaluator <- function(subset){
k<-5
set.seed(2)
ind <- sample(5, nrow(trainset), replace = TRUE)
results <- sapply(1:k, function(i){
train <- trainset[ind==i,]
test <- trainset[ind!=i,]
tree <- rpart(as.simple.formula(subset,"churn"), trainset)
error.rate <- sum(test$churn != predict(tree,test,type="class"))/nrow(test)
return(1- error.rate)
})
return(mean(results))
}
atrr.subset <- hill.climbing.search(names(trainset)[!names(trainset)%in% "churn"],evaluator)
f <- as.simple.formula(atrr.subset,"churn");f
churn ~ state + account_length + international_plan + voice_mail_plan +
number_vmail_messages + total_day_minutes + total_day_charge +
total_eve_calls + total_eve_charge + total_night_minutes +
total_night_calls + total_intl_calls + total_intl_charge +
number_customer_service_calls
<environment: 0x7f86b09bb4a0>
算法首先调用一个权重函数得到每个特征的权重值，权重评价指标是平均精确度下降importance.type = 1，除了上面用的随机森林，还可以使用chi.squared, information.gain。然后获取优化的特征子集，首先5折交叉验证评估特征子集的重要性，爬山搜索算法从原始特征集中选出优化的特征子集，也可以选择其他算法，比如 forward.search。还可以使用caret包进行特征筛选，据说这个包是个宝呀，包罗万象。
11.3 使用PCA进行降维
主成分分析是一种应用非常广泛的线性降维方法，适合数据集包含非常多的特征，并且特征间彼此冗余（相关的情况）。通过将特征集缩减成一小部分能代表原始特征集最主要变化的主要特征分量，实现高维数据到低维数据空间的映射。
# PCA
data("swiss")
swiss <- swiss[,-1]
# 两个参数分别是让变量向中心靠近，并使数据标准化
swiss.pca <- prcomp(swiss, center = TRUE, scale = TRUE)
swiss.pca
Standard deviations (1, .., p=5):
[1] 1.6228065 1.0354873 0.9033447 0.5592765 0.4067472
Rotation (n x k) = (5 x 5):
PC1 PC2 PC3 PC4 PC5
Agriculture 0.52396452 -0.25834215 0.003003672 -0.8090741 0.06411415
Examination -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
Education -0.49150243 0.19028476 0.539337412 -0.3321615 0.56656945
Catholic 0.38530580 0.36956307 0.725888143 0.1007965 -0.42176895
Infant.Mortality 0.09167606 0.87197641 -0.424976789 -0.2154928 0.06488642
summary(swiss.pca)
Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation 1.6228 1.0355 0.9033 0.55928 0.40675
Proportion of Variance 0.5267 0.2145 0.1632 0.06256 0.03309
Cumulative Proportion 0.5267 0.7411 0.9043 0.96691 1.00000
predict(swiss.pca, newdata=head(swiss,1))
PC1 PC2 PC3 PC4 PC5
Courtelary -0.9390479 0.8047122 -0.8118681 1.000307 0.4618643
特征选择过程中会去掉一些彼此关联但有价值的特征，需要在特征制取过程中考虑将这些特征综合到单特征中，PCA采用正交变换将彼此有关联的特征转化为主成分，以便我们确定方差趋势。算法主要包括以下步骤：1）找到平均向量的数据点；2）计算 协方差矩阵；3）计算特征向量；4）对特征向量排序并选择前k个特征向量；5）构建特征向量矩阵；最后，将数据样本转换成新的子集。拓展princomp是另一个高不成分分析函数，与上面的 prcomp采用奇异值分解不同，采用相关矩阵或协方差矩阵的特征值计算方法，一般更习惯用后者。
# princomp
swiss.princomp <- princomp(swiss, center = TRUE, scale = TRUE)
swiss.princomp
Call:
princomp(x = swiss, center = TRUE, scale = TRUE)
summary(swiss.princomp)
Standard deviations:
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
42.896335 21.201887 7.587978 3.687888 2.721105
5 variables and 47 observations.
Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Standard deviation 42.8963346 21.2018868 7.58797830 3.687888330
Proportion of Variance 0.7770024 0.1898152 0.02431275 0.005742983
Cumulative Proportion 0.7770024 0.9668177 0.99113042 0.996873399
Comp.5
Standard deviation 2.721104713
Proportion of Variance 0.003126601
Cumulative Proportion 1.000000000
predict(swiss.princomp,swiss[1,])
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
Courtelary -38.95923 -20.40504 -12.45808 -4.713234 1.46634
以上两个函数均来自stats包，还可以使用psych包中的principal函数进行：
# psych包principal
install.packages("psych")
install.packages("GPArotation")
library(psych)
swiss.principal <- principal(swiss, nfactors = 5, rotate = "none")
swiss.principal
Principal Components Analysis
Call: principal(r = swiss, nfactors = 5, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
PC1 PC2 PC3 PC4 PC5 h2 u2 com
Agriculture -0.85 -0.27 0.00 0.45 -0.03 1 0.00000000000000056 1.8
Examination 0.93 -0.01 -0.04 0.24 0.29 1 0.00000000000000278 1.3
Education 0.80 0.20 0.49 0.19 -0.23 1 0.00000000000000222 2.1
Catholic -0.63 0.38 0.66 -0.06 0.17 1 0.00000000000000111 2.8
Infant.Mortality -0.15 0.90 -0.38 0.12 -0.03 1 0.00000000000000067 1.5
PC1 PC2 PC3 PC4 PC5
SS loadings 2.63 1.07 0.82 0.31 0.17
Proportion Var 0.53 0.21 0.16 0.06 0.03
Cumulative Var 0.53 0.74 0.90 0.97 1.00
Proportion Explained 0.53 0.21 0.16 0.06 0.03
Cumulative Proportion 0.53 0.74 0.90 0.97 1.00
Mean item complexity = 1.9
Test of the hypothesis that 5 components are sufficient.
The root mean square of the residuals (RMSR) is 0
with the empirical chi square 0 with prob < NA
Fit based upon off diagonal values = 1
11.4 使用scree测试确定主成分数
Kaiser方法、scree(碎石测试）和依据挑选规则使用解释变量比例都可以。碎石测试的主要目的是将主成分结果以碎石图方式表达，从图中找到引起曲线斜率变化最快的因素。
# scree
screeplot(swiss.pca, type = "barplot")
screeplot(swiss.pca, type = "line")
主成分为2时，斜率变化最快。也可以使用nfactors以并行分析非图形方式作Cattell碎石来测试。
# nfactors
install.packages("nFactors")
library(nFactors)
ev <- eigen(cor(swiss))
ap <- parallel(subject=nrow(swiss), var=ncol(swiss), rep=100, cent=.05)
nS <- nScree(x=ev$values, aparallel = ap$eigen$qevpea)
jpeg(filename = 'test.jpg',width = 960,height = 480)
plot(nS)
dev.off()
碎石测试的非图形化实现
使用Kaider方法确定主成分数
# 标准差
swiss.pca$sdev
[1] 1.6228065 1.0354873 0.9033447 0.5592765 0.4067472
# 方差
swiss.pca$sdev^2
[1] 2.6335008 1.0722340 0.8160316 0.3127902 0.1654433
# 选择方差大于1的主成分
which(swiss.pca$sdev^2>1)
[1] 1 2
# 碎石图
screeplot(swiss.pca,type = "line")
abline(h=1, col="red", lty=3)
碎石图方差>1的主成分
11.6 主成分散点图可视化多元变量
biplot绘制数据与原始特征在前两个主成分上的投影图
# 散点图
plot(swiss.pca$x[,1], swiss.pca$x[,2], xlim = c(-4,4))
text(swiss.pca$x[,1], swiss.pca$x[,2], rownames(swiss.pca$x),
cex=0.7, pos = 4, col = "red")
biplot(swiss.pca)
biplot绘制数据及原始特征在前两个主成分上的投影，农业高，教育和检查低的省份在PC1上得分高；婴儿死亡率高，农业低的省份在主成分PC2上得分较高。
# ggbiplot
install.packages("devtools")
devtools::install_github("vqv/ggbiplot")
library(ggbiplot)
g <- ggbiplot(swiss.pca, obs.scale = 1, var.scale = 1,
eclispe = TRUE, circle = TRUE)
plot(g)
ggbiplot图
11.7 使用MDS进行降维分析
多维尺度分析通过图形方式展示多个对象之间的相似或相异程度距离），多维是指映射到一维、二维或多维空间表达CF全家人相对距离，一般使用一或二维空间。
# 计量MDS
swiss.dist <- dist(swiss)
swiss.mds <- cmdscale(swiss.dist,k=2)
plot(swiss.mds[,1], swiss.mds[,2], type = 'n', main = 'cmdscale(stats)')
text(swiss.mds[,1], swiss.mds[,2],rownames(swiss), cex=0.5, xpd=TRUE)
MDS二维散点图，好像没点？
分成计量和非计量两类，前者是主要考虑如何保证降维后各对象之间的距离尽可能接近它们在原始空间的距离，后者则假设两个空间中对象的距离排名已知，而且变换后排名不变。
# 非计算MDS
library(MASS)
swiss.nmds <- isoMDS(swiss.dist, k=2)
initial value 2.979731
iter 5 value 2.431486
iter 10 value 2.343353
final value 2.338839
converged
# 非计量MDS
library(MASS)
swiss.nmds <- isoMDS(swiss.dist, k=2)
plot(swiss.nmds$points, type = 'n', main = 'isoMDS(MASS)')
text(swiss.nmds$points, rownames(swiss), cex=0.5, xpd=TRUE)
isoMDS
# Shepard曲线
swiss.sh <- Shepard(swiss.dist, swiss.mds)
plot(swiss.sh, pch='.')
lines(swiss.sh$x, swiss.sh$yf, type = "S")
通常情况下选择二维欧氏空间作为目标空间，也可以选择其他距离或更高维空间。计量MDS方法使用基于计量的输入矩阵，也称主坐标分析，首先将距离转化为相似度；而非计量方法适合于顺序尺度的数据，仅考虑向量之间的距离排名，采用不可计量的输入矩阵。扩展将MDS对象作为一个图形对待是另一种可视化处理方法。
# MDS当图形
install.packages("igraph")
library(igraph)
swiss.sample <- swiss[1:10,]
g <- graph.full(nrow(swiss.sample))
V(g)$label <- rownames(swiss.sample)
layout <- layout.mds(g, dist = as.matrix(dist(swiss.sample)))
plot(g, layout=layout, vertex.size=3)
可以通过将投影维度绘制在一个散点图中比较MDS和PCA的差异，如果MDS采用欧氏距离，投影维度将与PCA完全一致。
# MDS Vs PCA
swiss.dist <- dist(swiss)
swiss.mds <- cmdscale(swiss.dist, k=2)
plot(swiss.mds[,1], swiss.mds[,2], type = 'n')
text(swiss.mds[,1], swiss.mds[,2], rownames(swiss), cex=0.9, xpd=TRUE)
swss.pca <- prcomp(swiss)
text(-swiss.pca$x[,1], -swiss.pca$x[,2], rownames(swiss),
col='blue', adj=c(0.2, -0.5), cex=0.9, xpd=TRUE)
MDS VS PCA，我的结果咋不一致呢？
11.8 使用SVD进行降维
奇异值分解是矩阵分解的一种形式，可以将一个矩阵分解为两个正交矩阵和一个对角矩阵，原始矩阵可由这三个矩阵相乘得到。可以帮助去掉那些从线性代数角度观察存在线性相关冗余的矩阵，可以应用在特征筛选，图像处理和聚类等。
# SVD
swiss.svd <- svd(swiss)
plot(swiss.svd$d^2/sum(swiss.svd$d^2), type = "l", xlab = "Singular vector",
ylab = "Variance explained")
plot(cumsum(swiss.svd$d^2/sum(swiss.svd$d^2)), type = "l", xlab = "Singular vector",
ylab = "Cumulative Variance explained")
# 只使用奇异向量重构数据
swiss.recon <- swiss.svd$u[,1] %*% diag(swiss.svd$d[1],
length(1), length(1)) %*% t(swiss.svd$v[,1])
par(mfrow=c(1,2))
image(as.matrix(swiss), main="Swiss data Image")
image(swiss.recon, main= "Reconstructed Image")
SVD是一类分解实数或复数矩阵的常见方法，PCA可以被看成SVD的一种特例：
svd.m <- svd(scale(swiss))
svd.m$v
[,1] [,2] [,3] [,4] [,5]
[1,] 0.52396452 -0.25834215 0.003003672 -0.8090741 0.06411415
[2,] -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
[3,] -0.49150243 0.19028476 0.539337412 -0.3321615 0.56656945
[4,] 0.38530580 0.36956307 0.725888143 0.1007965 -0.42176895
[5,] 0.09167606 0.87197641 -0.424976789 -0.2154928 0.06488642
pca.m <- prcomp(swiss, scale=TRUE)
pca.m$rotation
PC1 PC2 PC3 PC4 PC5
Agriculture 0.52396452 -0.25834215 0.003003672 -0.8090741 0.06411415
Examination -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
Education -0.49150243 0.19028476 0.539337412 -0.3321615 0.56656945
Catholic 0.38530580 0.36956307 0.725888143 0.1007965 -0.42176895
Infant.Mortality 0.09167606 0.87197641 -0.424976789 -0.2154928 0.06488642
两个矩阵基本相同。
11.9 使用SVD进行图像压缩
[图片上传失败...(image-be0ae8-1639570485003)] 图像压缩领域应用最为广泛的标准测试图像，花花公子当年的模特图呀！
# 图像压缩
download.file("lena512.bmp",url="https://gitee.com/zd200572/jpeg2bmp/raw/master/lenna512.bmp")
install.packages("bmp")
library(bmp)
lenna <- read.bmp("lena512.bmp")
lenna <- t(lenna[,,1])[,nrow(lenna[,,1]):1] # 导入时发生了旋转，转回去
image(lenna)
lenna.svd <- svd(scale(lenna))
plot(lenna.svd$d^2/sum(lenna.svd$d^2), type = "l", xlab = "Singular vector",
ylab = "Variance explained")
# 重构
length(lenna.svd$d)
[1] 512
min(which(cumsum(lenna.svd$d^2/sum(lenna.svd$d^2))>0.9))
[1] 17
lenna_compression <- function(dim) {
u<- as.matrix(lenna.svd$u[,1:dim])
v <- as.matrix(lenna.svd$v[,1:dim])
d <- as.matrix(diag(lenna.svd$d)[1:dim, 1:dim])
image(u%*%d%*%t(v))
}
lenna_compression(17)
min(which(cumsum(lenna.svd$d^2/sum(lenna.svd$d^2))>0.99))
# [1] 93
不知为啥，读什么图片都是负片呢？先继续：
11.10 使用ISOMAP进行非线性降维
ISOMAP属于流形学习方法，支持线性空间到非线性数据结构的转换，与MDS类似，它也能够以图形方式展现对象之间的相似性或相异性（距离），不过，由于数据采用非线性结构表示，以几何距离代替MDS中有欧氏距离。
# ISOMAP
install.packages("RnavGraphImageData")
library(RnavGraphImageData)
install.packages("vegan")
library(vegan)
data(digits)
sample.digit <- matrix(digits[,3000], ncol = 16, byrow=FALSE)
image(t(sample.digit)[,nrow(sample.digit):1])
# 随机选300个数字作为样本点
set.seed(2)
digit.idx <- sample(1:ncol(digits), size = 600)
digit.select <- digits[,digit.idx]
# vegdist计算对象之间的相异度
digits.Transpose <- t(digit.select)
digit.dist <- vegdist(digits.Transpose, method="euclidean")
# idomap降维
digit.isomap <- isomap(digit.dist, k=8, ndim = 6, fragmentedOK=TRUE) # 最后一个参数，最大连接可分析
plot(digit.isomap)
digit.st <- spantree(digit.dist)
digit.plot <- plot(digit.isomap, main="isomap k=8")
lines(digit.st, digit.plot, col = 'red')
ISOMAP是一种等距映射非线性降维方法，如果将计量MDS方法中数据点间成对的欧氏距离替换成邻接图间的测地距离，就可以将ISOMAP当做计量MDS方法的扩展。算法分为4步：确定邻近点，构建邻接图，计算最短路径和MDS分析找到数据间的低维嵌入。
扩展可以将RnavGraph包将图形作为数据浏览的基础方式来实现高维数据的可视化。
# 包的安装十分的不容易呀
install.packages("devtools")
BiocManager::install("graph")
devtools::install_git('https://gitee.com/zd200572/RnavGraph')
# bash git clone https://gitee.com/zd200572/RnavGraph.git
library(RnavGraph)
data(digits)
digit.group <- rep(c(1:9,0), each=1100)
digit.ng_data <- ng_data(name = "ISO_digits",
data = data.frame(digit.isomap$points),
shortnames = paste('i', 1:6, sep = ''),
group = digit.group[digit.idx],
labels = as.character(digit.group[digit.idx]))
V <-shortnames(digit.ng_data)
G <- completegraph(V)
LG <- linegraph(G)
LGnot <- complement(LG)
ng.LG <- ng_graph(name = "3D Transition", graph = LG)
ng.LGnot <- ng_graph(name = "4D Transition", graph = LGnot)
ng.i.digits = ng_image_array_gray("UPS Handwriting Digits",
digit.select,16,16,invert = TRUE,
img_in_row = FALSE)
vizDigits1 <- ng_2d(data = digit.ng_data, graph = ng.LG, images = ng.i.digits)
vizDigits2 <- ng_2d(data = digit.ng_data, graph = ng.LGnot, images = ng.i.digits)
nav <- navGraph(data = digit.ng_data, graph = list(ng.LG, ng.LGnot),
viz = list(vizDigits1, vizDigits2))
11.11 使用局部线性嵌入法进行非线性降维
LLE算法是PCA算法的扩展，通过嵌入高维空间内的流形映射到低维空间来实现数据压缩。ISOMAP是全局性非线性降维，LLE主要是局部母性降维算法，假设每个数据点可以由k个邻近点的母性组合构成，映射后能保持原来的数据性质。
# LLE
install.packages("lle")
library(lle)
data("lle_scurve_data")
X <- lle_scurve_data
results <- lle(X=X, m=2, k=12,id=TRUE) # 保留维度为2，邻近点设为12，参数id TRUE
finding neighbours
calculating weights
intrinsic dim: mean=2.47875, mode=2
computing coordinates
str(results)
List of 4
$ Y : num [1:800, 1:2] -1.586 -0.415 0.896 0.513 1.477 ... #嵌入数据
$ X : num [1:800, 1:3] 0.955 -0.66 -0.983 0.954 0.958 ... # 输入数据
$ choise: NULL # 保留数据的索引向量？
$ id : num [1:800] 3 3 2 3 2 2 2 3 3 3 ... # 输入数据的维度
plot(results$Y, main = "embeded data", xlab = expression(y[1]),
ylab = expression(y[2]))
plot_lle(results$Y, X, FALSE, col = 'red', inter = TRUE) # 最后一个参数是交互式
LLE是一种非线性降维算法，基于它我们可以得到高维数据在低维空间保持原有数据邻近嵌入关系的映射。算法主要分成三步：计算每个点的k个邻近，然后计算每个邻近点的权值，使得每个点都能最优地由其邻近点组合重构，即残差和最小。
扩展还可以选择RDRTollbox包实现非线性降维，支持ISOMAP和LLE算法。
BiocManager::install("RDRToolbox")
library(RDRToolbox)
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2021-12-28，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读