当前位置：首页 > article >正文

PCA降维算法

article 2024/10/7 10:52:08

假设原数据 $X$ 是由·m·个·n·维向量组成，现在要将数据从n维降到k维：

一、如何计算：

流程如下：

0.定义数据：
1.中心化：将 X 的每一维进行零均值化，即每个对应维分量都减去这一维度的均值
2.求出协方差矩阵 $C=\frac{1}{m}X^TX$ ：
3.求出协方差矩阵的特征值及对应的特征向量：
4.将特征向量按对应特征值大小从左到右按列排列成矩阵，取前 k 列组成矩阵 P
5. $Y = XP$ 即为降维到 k 维后的数据

0.定义数据：

如：现在有3个4维向量，我们要将每个向量降到2维：（3 * 4）->（3 * 2）

data = np.array([[1,5,3,1],
                 [2,3,4,4],
                 [6,7,8,1]])
print("data.shape:",data.shape)

输出：

data.shape: (3, 4)

1.中心化：

将 X 的每一维进行零均值化，即每个对应维分量都减去这一维度的均值

对于第一维：μ=(1+2+6)/3=3 ，均值为3，更新第一维的分量为-2，-1，3
对于第二维：μ=(5+3+7)/3=5 ，均值为5，更新第一维的分量为0，-2，2
对于第三维：μ=(3+4+8)/3=5 ，均值为5，更新第一维的分量为-2，-1，3
对于第四维：μ=(1+4+1)/3=2 ，均值为2，更新第一维的分量为-1，2，-1

# 中心化
Mean=data.mean(axis=0)

print("data:\n",data)
print("Mean:\n",Mean)
data = data - Mean
print("newdata:\n",data)

输出：

data:
[[1 5 3 1]
[2 3 4 4]
[6 7 8 1]]
Mean:
[3. 5. 5. 2.]
newdata:
[[-2. 0. -2. -1.]
[-1. -2. -1. 2.]
[ 3. 2. 3. -1.]]

2.求出协方差矩阵 $C=\frac{1}{m}X^TX$ ：

# 注意：计算协方差矩阵,求出来应该是4*4的矩阵，因为有4个特征
#打印维数
cov = np.dot(data.T, data) / (m) 
print("cov:\n",cov)

输出：

cov:
[[ 4.66666667 2.66666667 4.66666667 -1. ]
[ 2.66666667 2.66666667 2.66666667 -2. ]
[ 4.66666667 2.66666667 4.66666667 -1. ]
[-1. -2. -1. 2. ]]

3.求出协方差矩阵的特征值及对应的特征向量：

# 求出协方差矩阵的特征值及对应的特征向量
eig_val, eig_vec = np.linalg.eig(cov)
print("特征值:\n",eig_val)
print("特征向量:\n",eig_vec)

特征值:
[1.15092498e+01 2.49075025e+00 9.25112518e-17 8.09374965e-16]
特征向量:
[[-0.62111719 0.30159181 0.07152408 -0.42472814]
[-0.42435516 -0.48839262 0.72580926 0.69772981]
[-0.62111719 0.30159181 -0.36184778 0.14563622]
[ 0.21988535 0.76128668 0.58064741 0.55818385]]

每个特征值的特征向量对应下面矩阵中的一列，而不是一行。

4.将特征向量按对应特征值大小从左到右按列排列成矩阵，取前 k 列组成矩阵 P

# 将特征向量按对应特征值大小从左到右按列排列成矩阵，取前 k 列组成矩阵 P
k=2
P = eig_vec[:, np.argsort(-eig_val)[:k]]
print("基向量矩阵:\n",P)

基向量矩阵:
[[-0.62111719 0.30159181]
[-0.42435516 -0.48839262]
[-0.62111719 0.30159181]
[ 0.21988535 0.76128668]]

注意：
[-0.62111719， -0.42435516 ，-0.62111719 ， 0.21988535]对应一个特征向量，即是按列看的。

5. $Y = XP$ 即为降维到 k 维后的数据

降维后的结果:

# Y = XP
Y = np.dot(data, P)
print("降维后的结果:\n",Y)

输出：

Y= [[ 2.26458342 -1.96765391]
[ 2.5307154 1.89617499]
[-4.79529882 0.07147892]]

两个矩阵相乘的意义是将左边矩阵中的每一行向量变换到右边矩阵中以每一列列向量为基所表示的空间中去（这里左边矩阵的每一行代表一个数据点，右边矩阵的每一列代表一个基）

$Y = XP$ ，即将X里面的点投影到P所构成的空间中去

与PCA函数的计算结果进行对比：结果一致

# # PCA 降维到 2 维
pca = PCA(n_components=2)
reduced_data = pca.fit_transform(data)

print(reduced_data)

[[-2.26458342 1.96765391]
[-2.5307154 -1.89617499]
[ 4.79529882 -0.07147892]]

二、判断需要降到多少维：

通过计算每个PC的方差占总方差的比例，来确定需要降到多少维，既可以减少数据的存储开销，又能尽量减少数据损失。

PCA算法的实现可以分为最大可分性和最近重构性，前者更容易编程实现，显然我们使用的也是前者。

目标在于使投影后的方差最大化，但是经过数学推导后可以发现PC对应的方差其实就是协方差矩阵的特征值，所以我们可以通过协方差矩阵的特征值占比来确定需要降到多少维。


# -------------判断需要降到几维----------------

# 计算每个主成分的方差贡献率，即特征值的大小与特征值之和的比值
# 方差贡献率百分比，保留两位小数并添加百分号
var_ratio = eig_val / eig_val.sum()
var_ratio_percent = ["{:.2f}%".format(vr * 100) for vr in var_ratio]
print("方差贡献率百分比:\n", var_ratio_percent)


# 绘制方差贡献率的条形图
plt.figure(figsize=(10, 6))
plt.bar(range(len(var_ratio)), [vr * 100 for vr in var_ratio], tick_label=['PC1', 'PC2', 'PC3', 'PC4'])
plt.xlabel('λ')
plt.ylabel('rate (%)')
plt.title('rate of λ')
plt.show()

画图如下：
在这里插入图片描述
从图中也可以看出，PC1和PC2几乎占据了所有的信息，所以降到2维就是最好的选择。

由于计算机在计算浮点数时会产生一些误差，所以出现一些奇怪的值也不必惊慌。

三、降维后是否还能恢复：

使用了 PCA（主成分分析）降维算法后，一般无法完全恢复原始数据。这是因为在降维的过程中，我们丢弃了某些主成分（通常是那些对应较小特征值的主成分），从而丢失了一部分信息。

但是可进行近似重构，使用保留的主成分对数据进行重构，得到对原始数据的近似，但会有信息损失和误差。

代码

# -*- coding: utf-8 -*-
import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# 示例数据矩阵
data = np.array([[1,5,3,1],
                 [2,3,4,4],
                 [6,7,8,1]])

m=3 # 样本数
n=4 # 特征数

print("data.shape:",data.shape)

# 去中心化
Mean=data.mean(axis=0)

print("data:\n",data)
print("Mean:\n",Mean)
data = data - Mean
print("newdata:\n",data)

# 计算协方差矩阵,求出来应该是4*4的矩阵，因为有4个特征
#打印维数
cov = np.dot(data.T, data) / (m) 
print("cov:\n",cov)

# 求出协方差矩阵的特征值及对应的特征向量
eig_val, eig_vec = np.linalg.eig(cov)
print("特征值:\n",eig_val)
print("特征向量:\n",eig_vec)


# 将特征向量按对应特征值大小从上到下按行排列成矩阵，取前 k 行组成矩阵 P
k=2
P = eig_vec[:, np.argsort(-eig_val)[:k]]
print("基向量矩阵:\n",P)


# Y = XP
Y = np.dot(data, P)
print("降维后的结果:\n",Y)


# # PCA 降维到 2 维
pca = PCA(n_components=2)
reduced_data = pca.fit_transform(data)

print(reduced_data)

# -------------判断需要降到几维----------------

# 计算每个主成分的方差贡献率，即特征值的大小与特征值之和的比值
# 方差贡献率百分比，保留两位小数并添加百分号
var_ratio = eig_val / eig_val.sum()
var_ratio_percent = ["{:.2f}%".format(vr * 100) for vr in var_ratio]
print("方差贡献率百分比:\n", var_ratio_percent)


# 绘制方差贡献率的条形图
plt.figure(figsize=(10, 6))
plt.bar(range(len(var_ratio)), [vr * 100 for vr in var_ratio], tick_label=['PC1', 'PC2', 'PC3', 'PC4'])
plt.xlabel('λ')
plt.ylabel('rate (%)')
plt.title('rate of λ')
plt.show()