Algorithm[1] PCA主成分分析原理剖析与实现

悬铃木下的青春

已于 2022-09-28 21:16:50 修改

阅读量419

点赞数

分类专栏：算法笔记文章标签：机器学习算法人工智能

于 2022-09-17 21:18:29 首次发布

本文链接：https://blog.csdn.net/m0_46500149/article/details/126898677

版权

算法笔记专栏收录该内容

1 篇文章 0 订阅

订阅专栏

定义

通过正交变换将一组可能存在相关性的变量转换为一组线性无关的变量，转换后的这组变量叫主成分。

本质

在n维向量空间中，寻找使原样本方差最大的基向量组，并将数据投影到新的向量空间中。
我们无论是在处理、分析或是存储数据时，都希望用相对少的特征去表征样本的性态。
一方面，我们希望不同特征之间都是相互独立的，这样可以减少数据的冗余，但特征之间总会呈现出一定的相关性，当两个特征强相关时，我们事实上可以用一个新的特征去代替这两个特征，从而实现降维。
另一方面，方差越大，说明样本之间的差异更加明显，其蕴含的信息也就越多，价值也相应越大。这决定了我们确定基向量的策略，这也正是PCA的思想所在。

主要作用

对于多特征问题，可以有效降维，实现特征筛选。
消除数据冗余，去除部分噪声，降低处理难度。
将高维数据转为二维或三维问题，便于数据可视化。
构造筛选回归模型因子，简化回归模型。

数学推导

Q1:协方差矩阵计算公式

在这里插入图片描述

Q2:SVD降维原理

在这里插入图片描述

Q3:贡献率公式

$=\frac{\lambda_i}{\sum\lambda_i}$

算法步骤

已知样本 $X_{origin}=\{ x_1,x_2, ..,x_n \}$ ，组成 $m * n$ 维样本矩阵（即 $m$ 个样本，每个样本有 $n$ 个特征）

注意：由于数据可能存在缺失、异常或量纲不一致问题，这里默认数据已经预处理

特征分解法

对样本 $X_{origin}$ 去均值，得新样本矩阵 $X$
求协方差矩阵 $\frac{1}{m-1}X^TX$
协方差矩阵对角化
解得特征值对角阵 $\Lambda_{n*n}$ 、特征向量矩阵 $P_{n*n}$
计算各特征值贡献率，确定主成分个数 $k$
获得降维后的主成分分数矩阵 $Y_{m*k}=X_{m*n}P_{n*k}$

奇异值分解法（SVD）

散度矩阵 $X^TX$ 对角化（列降维）（有些算法可以跳过这一步）
解得奇异值矩阵 $\Sigma$ 和对应右奇异矩阵 $V$
计算奇异值贡献率，确定主成分个数 $k$
获得降维后的主成分分数矩阵 $Y_{m*k}=X_{m*n}V^{T}_{n*k}$

补充：若需要对行降维，则对散度矩阵 $XX^T$ 对角化获得左奇异矩阵，最终由 $Y_{k*n}=U_{k*m}X_{m*n}$ 得到降维后的矩阵

算法实现

Python

法一：自己撸脚本

注意：该函数有两种赋值方式 —— 除了样本矩阵X为必填参数，形参 percent 和 k 可以任选其一，不可以同时赋值，可能会发生溢出错误！同时缺省则默认输出1个主成分。

import numpy as np
'''
--------------
实现PCA算法函数
--------------
X为样本矩阵
percent为指定贡献率（默认值0）
k为主成分个数（默认值1）
'''
def pca(X, percent = 0, k = 1):
    #返回对应的样本数（行）和特征数（列）
    n_samples, n_features = X.shape
    #列表解析用法，循环计算每一列的样本均值，返回列表，再生成数组
    mean = np.array([np.mean(X[:,i]) for i in range(n_features)])
    #去均值
    norm_X = X - mean
    #计算协方差矩阵
    cov_matrix=np.dot(np.transpose(norm_X),norm_X)/(n_samples-1)
    #计算特征值（未按次序，复数域）和特征向量
    eig_val, eig_vec = np.linalg.eig(cov_matrix)
    #将特征值与特征向量一一对应
    eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(n_features)]
    #基于特征值降序排列（True为降序）
    eig_pairs.sort(reverse=True)
    #初始化累计贡献值
    Con = np.sum(eig_val[:k])/np.sum(eig_val)
    #计算累计贡献值达设定值需要多少个主成分
    while Con < percent:   
        k += 1 
        Con = np.sum(eig_val[:k])/np.sum(eig_val)
        
    #选择前k个特征向量
    feature = np.array([ele[1] for ele in eig_pairs[:k]])
    #feature变成了k*n的矩阵,须转置
    data = np.dot(norm_X,np.transpose(feature))
    return data

X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
print(pca(X,k=1)) #指定主成分个数解出主成分分数矩阵
print(pca(X,0.9)) #根据累计贡献率解出主成分分数矩阵
print(pca(X)) #全缺省，默认输出一个主成分

法二：调库

sklearn库实现了PCA算法，这里大致介绍一下这个函数用法：

sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)

n_components：所要保留的主成分个数（1,2,…（int））
copy:保存原始数据，默认保存
whiten:白化，使得每个特征具有相同的方差，默认False

通过创建一个PCA对象实例，我们可以对它进行以下操作：

方法：
- fit(X,y=None) X为输入数据，训练pca模型。无监督学习。
- transform(X) 将数据X转换成降维后的数据
- fit_transform(X) 用X来训练PCA模型，同时返回降维后的数据
- inverse_transform(new_X) 将降维后的数据转换成原始数据
属性：
- components_ ：返回特征向量矩阵，注意是按行排列。
- explained_variance_：返回特征值
- explained_variance_ratio_：返回所保留的n个成分各自的方差百分比。
- n_components_：返回所保留的成分个数n。

from sklearn.decomposition import PCA
'''
--------------------------------------------------------------------
sklearn包中decomposition模块实现了PCA算法
--------------------------------------------------------------------
'''
X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca = PCA(n_components=1)  #返回一个pca对象
pca_data = pca.fit_transform(X)  #训练
print(pca_data)

Matlab

法一：自己撸脚本

Matlab对于矩阵运算比较方便

%matlab实现pca
A = [-1 1;-2 -1;-3 -2;1 1;2 1;3 2];

%数据标准化处理
[a, b] = size(A,1,2); % a-样本数, b-特征数

%提高运算速度，给矩阵预分配内存(针对随迭代动态改变的)
SA = zeros(a,b);
DS = zeros(b,3);

%去均值
for i=1:b
    SA(:,i) = A(:,i)-mean(A(:,i));
end

% 计算协方差矩阵
CM = cov(SA);
%计算特征值和特征向量。
[V,D]=eig(CM); 

for j=1:b
    DS(j,1)=D(b+1-j,b+1-j);%将特征值按降序排序，放在DS第1列
end

for i=1:b
    DS(i,2)=DS(i,1)/sum(DS(:,1)); %计算贡献率，放在DS第2列
    DS(i,3)=sum(DS(1:i,1))/sum(DS(:,1));%计算累积贡献率，放在DS第3列
end

T=0.9;%主成分信息保留率

for K=1:b
    if DS(K,3)>=T
        Com_num=K; %达到90%即可，记录主成分个数
        break;
    end
end

%提高运算速度，给矩阵预分配内存(针对随迭代动态改变的)
PV = zeros(b,Com_num);
total_score = zeros(a,2);

%提取主成分对应的特征向量
for j=1:Com_num
    PV(:,j)=V(:,b+1-j); %V是对应的特征向量，按列排序
end

%计算各评价对象的主成分得分
new_score = SA*PV; %SA：15*8，PV：8*3 生成PCA降维后的new_score矩阵


%输出模型及结果报告
disp('特征值及其贡献率、累计贡献率：')
DS %贡献率表
disp('信息保留率T对应的主成分数与特征向量：')
Com_num %主成分个数
PV %对应的特征向量
disp('生成PCA降维后的主成分分数矩阵：')
new_score %分数矩阵

法二：调用函数

Matlab实现了PCA算法，这里大致介绍一下这个函数用法：

coeff = pca(X)
[coeff,score,latent] = pca(X)
[coeff,score,latent,~,explained] = pca(X); %不想输出的项可以用~

PCA 默认数据中心化，并使用奇异值分解 (SVD) 算法。
输入参数
NumComponents：可以指定主成分个数
输出参数
在 coeff 中返回大小为 $n\times n$ 的主成分系数矩阵，即特征向量矩阵，其中列按方差的降序排列。
在 score 中返回主成分分数，矩阵大小为 $m\times k$ 。
在 latent 中返回主成分方差，即协方差矩阵的特征值，以列向量形式返回。
在 explained 中返回特征值贡献率。

% matlab实现pca
X = [-1 1;-2 -1;-3 -2;1 1;2 1;3 2];
[coeff,score,latent]=pca(X)

总结

自己写的代码是通过求解协方差矩阵来实现PCA的，而Python的sklearn包以及Matlab内置函数中默认采用SVD策略实现PCA。
特征分解只能对一个维度进行降维，而SVD可以通过右奇异矩阵按列降维，亦可以通过左奇异矩阵按行降维。即降低冗余特征的同时还可以降低冗余样本。
SVD相比特征分解更加稳定，时间、空间开销更小，适合稀疏矩阵的降维。特征分解法需要计算协方差矩阵，且需要进行去均值处理，会破坏矩阵的稀疏性，增大存储负担。
SVD可以作为PCA算法中的一个辅助算法，亦可作为图像压缩、自然语言处理等领域的算法。