首页 > 分享 > PCA主成分分析去噪与降维

PCA主成分分析去噪与降维

花匠小妙招
2025-02-08 10:15

PCA主成分分析去噪与降维

（项目学习笔记及代码）

一、目的

1、降低维数（主要目的）

主成分分析 (Principal Component Analysis，简称 PCA)是最常用的一种降维方法。（具体可见西瓜书）

2、去噪

PCA也可以用作噪音的过滤，因为任何一个成分的变化影响都远远大于随机噪声的影响的，所以针对于噪声，各种的成分相对不受影响。即可以用主成分来重构原始数据。

3、原理假设 R d mathbb{R}^{d} Rd空间中有 m m m个样本点 : X = { x 1 , x 2 , . . . , x m } :X={x_1,x_2,...,x_m} :X={x1,x2,...,xm}。这些点即可以理解为我们数据的特征，如图片的像素，每一个数据 x i x_i xi都有 d d d维， x i ∈ R d x_i in mathbb{R}^{d} xi∈Rd。
x i = [ x i 1 x i 2 ⋅ ⋅ ⋅ x i d ] d × 1 (1) begin{aligned} x_i={left[ begin{matrix} x_i^1 \ x_i^2 \ cdot \ cdot \ cdot \ x_i^d \ end{matrix} right] }_{d times 1} tag{1} end{aligned}

begin{aligned} x_i={left[ begin{matrix} x_i^1 \ x_i^2 \ cdot \ cdot \ cdot \ x_i^d \ end{matrix} right] }_{d times 1} tag{1} end{aligned}

xi=⎣⎢⎢⎢⎢⎢⎢⎡xi1xi2⋅⋅⋅xid⎦⎥⎥⎥⎥⎥⎥⎤d×1(1)
于是有：
X = ( x 1 , x 2 , . . . , x m ) d × m T = [ x 1 1 x 2 1 ⋅ ⋅ ⋅ x 1 d x 1 2 x 2 2 ⋅ ⋅ ⋅ x 2 d ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ x m 1 x m 2 ⋅ ⋅ ⋅ x m d ] m × d (2) begin{aligned} X=(x_1,x_2,...,x_m)_{d times m}^mathrm{T}={left[ begin{matrix} x_1^1 & x_2^1 & cdot & cdot & cdot & x_1^d \ x_1^2 & x_2^2 & cdot & cdot & cdot & x_2^d \ cdot & cdot & cdot & & & cdot \ cdot & cdot & & cdot & & cdot \ cdot & cdot & & & cdot & cdot \ x_m^1 & x_m^2 & cdot & cdot & cdot & x_m^d \ end{matrix} right] }_{m times d} tag{2} end{aligned}

begin{aligned} X=(x_1,x_2,...,x_m)_{d times m}^mathrm{T}={left[ begin{matrix} x_1^1 & x_2^1 & cdot & cdot & cdot & x_1^d \ x_1^2 & x_2^2 & cdot & cdot & cdot & x_2^d \ cdot & cdot & cdot & & & cdot \ cdot & cdot & & cdot & & cdot \ cdot & cdot & & & cdot & cdot \ x_m^1 & x_m^2 & cdot & cdot & cdot & x_m^d \ end{matrix} right] }_{m times d} tag{2} end{aligned}

X=(x1,x2,...,xm)d×mT=⎣⎢⎢⎢⎢⎢⎢⎡x11x12⋅⋅⋅xm1x21x22⋅⋅⋅xm2⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅x1dx2d⋅⋅⋅xmd⎦⎥⎥⎥⎥⎥⎥⎤m×d(2)
其中 x i j x_i^j xij表示是第 i i i个样本的第 j j j维的特征，单独的列向量 x i x_i xi即是 d d d维的。每一行代表一个数据样本。首先对 x i x_i xi数据进行归一化（保证其均值为0）。
x ˉ = 1 m ∑ i = 1 m x i = [ x ˉ 1 x ˉ 2 ⋅ ⋅ ⋅ x ˉ d ] d × 1 (3) bar{x}= frac{1}{m} sum_{i=1}^m x_i={left[ ˉx1ˉx2⋅⋅⋅ˉxd right] }_{d times 1} tag{3} xˉ=m1i=1∑mxi=⎣⎢⎢⎢⎢⎢⎢⎡xˉ1xˉ2⋅⋅⋅xˉd⎦⎥⎥⎥⎥⎥⎥⎤d×1(3)
对所有的数据样本点列向量 x i x_i xi求对应元素的均值，得到d维的均值列向量 x ˉ bar{x} xˉ，然后用原始样本数据减去均值即得到了归一化后的结果。

x i = x i − x ˉ = [ x i 1 − x ˉ 1 x i 2 − x ˉ 2 ⋅ ⋅ ⋅ x i d − x ˉ d ] d × 1 (5) begin{aligned} x_i & =x_i-bar{x}={left[ begin{matrix} x_i^1 - bar{x}^1 \ x_i^2 - bar{x}^2 \ cdot \ cdot \ cdot \ x_i^d - bar{x}^d \ end{matrix} right] }_{d times 1} tag{5} end{aligned}

xi=xi−xˉ=⎣⎢⎢⎢⎢⎢⎢⎡xi1−xˉ1xi2−xˉ2⋅⋅⋅xid−xˉd⎦⎥⎥⎥⎥⎥⎥⎤d×1(5)

对数据样本 x i ∈ R d x_iin mathbb{R}^{d} xi∈Rd进行降维压缩得到结果 z i ∈ R d ′ z_iin mathbb{R}^{d'} zi∈Rd′， ( d ′ < d ) (d'<d) (d′<d)。由 x i x_i xi得到 z i z_i zi的过程可以如下表示，权重 W = { w 1 , w 2 , . . . , w d ′ } mathbf{W}={w_1,w_2,...,w_{d'} } W={w1,w2,...,wd′}作用下的矩阵线性变换。这就是我们降维后的结果，也是PCA工作的目标：
z i j = w j T ⋅ x i = [ w j 1 w j 2 ⋅ ⋅ ⋅ w j d ] d × 1 T ⋅ [ x 1 x 2 ⋅ ⋅ ⋅ x d ] d × 1 (6) z_i^j=w_j^mathrm{T} cdot x_i={left[ w1jw2j⋅⋅⋅wdj right] }_{d times 1}^mathrm{T} cdot {left[ x1x2⋅⋅⋅xd right] }_{d times 1} tag{6} zij=wjT⋅xi=⎣⎢⎢⎢⎢⎢⎢⎡wj1wj2⋅⋅⋅wjd⎦⎥⎥⎥⎥⎥⎥⎤d×1T⋅⎣⎢⎢⎢⎢⎢⎢⎡x1x2⋅⋅⋅xd⎦⎥⎥⎥⎥⎥⎥⎤d×1(6)
其中 z i j z_i^j zij是降低维度后的原始样本数据 x i x_i xi对应的第 j j j维的元素值；而 w j w_j wj是 d d d维的列向量。从最大可分（最大投影方差）的角度出发，同时当 ∥ w j ∥ 2 2 = 1 |w_j |_2^2=1 ∥wj∥22=1时， w j T ⋅ x i w_j^mathrm{T} cdot x_i wjT⋅xi即是样本数据点在该矩阵方向上得投影。因为 α alpha α到 β beta β的投影即是 ∣ α ∣ cos ⁡ θ |alpha|cos theta ∣α∣cosθ，而向量点乘有 α ⋅ β = ∣ α ∣ ∣ β ∣ cos ⁡ θ alpha cdot beta=|alpha||beta|cos theta α⋅β=∣α∣∣β∣cosθ，当 ∣ β ∣ = 1 |beta|=1 ∣β∣=1时，点乘就是投影了。我们的目标就是找到恰当的 W mathbf{W} W将这个投影 w j T ⋅ x i w_j^mathrm{T} cdot x_i wjT⋅xi的方差最大化，也就是让样本点在新的基准下划分得越分开越好。于是有样本方差来度量：
（均值已经归一化为0了）
S = 1 m − 1 ∑ i = 1 m ( w j T x i − w j T x ˉ ) 2 = 1 m − 1 ∑ i = 1 m ( w j T ( x i − x ˉ ) ) 2 = 1 m − 1 ∑ i = 1 m ( w j T ( x i − x ˉ ) ) ⋅ ( w j T ( x i − x ˉ ) ) T = 1 m − 1 ∑ i = 1 m w j T ( x i − x ˉ ) ⋅ ( x i − x ˉ ) T w j → w j T ( ∑ i = 1 m ( x i − x ˉ ) ⋅ ( x i − x ˉ ) T ) w j → w j T ( Cov ( x ) ) w j → w j T S w j (7) begin{aligned} S &= frac{1}{m-1} sum_{i=1}^m (w_j^mathrm{T}x_i-w_j^mathrm{T}bar{x})^2\ &= frac{1}{m-1} sum_{i=1}^m (w_j^mathrm{T}(x_i-bar{x}))^2\ &= frac{1}{m-1} sum_{i=1}^m (w_j^mathrm{T}(x_i-bar{x})) cdot (w_j^mathrm{T}(x_i-bar{x}))^mathrm{T} \ &= frac{1}{m-1} sum_{i=1}^m w_j^mathrm{T}(x_i-bar{x}) cdot (x_i-bar{x})^mathrm{T}w_j \ & rightarrow w_j^mathrm{T} left(sum_{i=1}^m (x_i-bar{x}) cdot (x_i-bar{x})^mathrm{T} right)w_j \ & rightarrow w_j^mathrm{T} left(text{Cov} (x) right)w_j \ & rightarrow w_j^mathrm{T}S w_j \ tag{7} end{aligned} S=m−11i=1∑m(wjTxi−wjTxˉ)2=m−11i=1∑m(wjT(xi−xˉ))2=m−11i=1∑m(wjT(xi−xˉ))⋅(wjT(xi−xˉ))T=m−11i=1∑mwjT(xi−xˉ)⋅(xi−xˉ)Twj→wjT(i=1∑m(xi−xˉ)⋅(xi−xˉ)T)wj→wjT(Cov(x))wj→wjTSwj(7)
所以投影的方差最后转换成了样本矩阵的协方差。于是得到优化目标：
arg max ⁡ w j w j T S w j s.t. w j ⋅ w j T = 1 (8) argmax_{w_j} w_j^mathrm{T}S w_j \ text{s.t. } w_j cdot w_j^mathrm{T}=1 tag{8} wjargmaxwjTSwjs.t. wj⋅wjT=1(8)
利用拉格朗日数乘法：
L ( w j , λ ) = w j T S w j + λ ( w j ⋅ w j T − 1 ) (9) L(w_j,lambda)=w_j^mathrm{T}S w_j + lambda (w_j cdot w_j^mathrm{T}-1) tag{9} L(wj,λ)=wjTSwj+λ(wj⋅wjT−1)(9)
再对每一个 w j w_j wj求偏微分 ∂ L ∂ w j = 0 frac{partial L}{partial w_j}=0 ∂wj∂L=0：
S w j − λ w j = 0 S w j = λ w j (10) Sw_j-lambda w_j=0\ Sw_j=lambda w_j \ tag{10} Swj−λwj=0Swj=λwj(10)
于是乎，形如此形式有： w j w_j wj即是特征向量，而 λ lambda λ是特征值。
再带回优化目标有：
w j T S w j = w j T λ w j = λ (11) w_j^mathrm{T}S w_j=w_j^mathrm{T}lambda w_j = lambda tag{11} wjTSwj=wjTλwj=λ(11)
这意味这最大化的结果值（方差最大化的结果值）就是数据样本协方差矩阵的特征值。 可以理解到特征值在某种意义上表达的是方差的意义。
依次类推，此后得到的每一个的 w j w_j wj还需要保证相互之间点乘为0是正交的。于是一般即可以理解为找到数据样本 X X X的协方差矩阵 S S S，是一个实对称矩阵，那么其一定是可以相似对角化的，找到其最大的特征值 λ lambda λ即是最大的投影方差，其对应的特征向量 w j w_j wj即是线性乘法因子。而最后选择最大的特征值即可。一般来说目前都直接奇异值分解法求解，因为奇异值分解数据矩阵 X X X的过程也是会构造出 X X T XX^mathrm{T} XXT，也就是协方差矩阵。

二、代码

1、本文代码是通过加载原始的matlab 中的一个数组文件作为数据集的（数据集可以是任意的数据，但注意最后用plt库绘图时，注意下标的匹配）
2、其中有些地方代码已经注释的很清楚了。
3、若进行降噪即最后重构（restore）即可，若进行降维，那么经过转换后（transform）的即是新的维度下的结果。
4、代码引用了部分搜索到的，不过有些确实找不到了而没有引用！

1、sklearn库方式

# -*- coding: utf-8 -*- """ Created on Tue Jan 8 10:39:09 2019 @author: QSY """ import matplotlib.pyplot as plt import scipy.io as sci from sklearn.decomposition import PCA import numpy as np """Read data from *.mat including two kinds""" data_tarin = sci.loadmat('data_train.mat') # data_test = sci.loadmat('data_test.mat') data_train_noisy = data_tarin['noiseSignal'] data_train_theory = data_tarin['clearSignal'] #选取数据进行操作 X=data_train_noisy[0:1,:] #进行行列转换 X=list(map(list, zip(*X))) #创建横坐标 p=[] for i in range(434): j=i p.append(j) def ReadReductionofcomponents(X): pca=PCA(n_components=1) pca.fit(X) #进行操作转换后的结果 x_reduction =pca.transform(X) #进行复原恢复后的结果 x_restore=pca.inverse_transform(x_reduction) x_restore=x_restore[:,0:1] #X仅选择某一列的数据进行操作 X=np.array(X) X=X[:,0:1] #进行线绘制 plt.plot(p, x_restore, color='r', linewidth=0.3, alpha=0.6,label="PCA") plt.plot(p, X, color='b', linewidth=0.2, alpha=0.6,label="init") plt.legend()#显示图例名字 #pdf保存 plt.savefig('1.png') plt.show() ReadReductionofcomponents(X)

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455 2、数学方式

# -*- coding: utf-8 -*- """ Created on Tue Jan 8 16:24:57 2019 @author: QSY """ import numpy as np import scipy.io as sci import matplotlib.pyplot as plt """Read data from *.mat including two kinds""" data_tarin = sci.loadmat('data_train.mat') other_data=sci.loadmat('pca.mat') # data_test = sci.loadmat('data_test.mat') data_train_noisy = data_tarin['noiseSignal'] data_train_theory = data_tarin['clearSignal'] perfect_data=other_data['U'] #选取数据进行操作 X=data_train_noisy[0:3,:] X=list(map(list, zip(*X))) X=perfect_data #零均值化 def zeroMean(dataMat): #按列求均值，即求各个特征的均值 meanVal=np.mean(dataMat,axis=0) newData=dataMat-meanVal return newData,meanVal def pca(dataMat,n): newData,meanVal=zeroMean(dataMat) covMat=np.cov(newData,rowvar=0) #求协方差矩阵 eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特征值和特征向量,特征向量是按列放的，即一列代表一个特征向量 eigValIndice=np.argsort(eigVals) #对特征值从小到大排序 n_eigValIndice=eigValIndice[-1:-(n+1):-1] #最大的n个特征值的下标 n_eigVect=eigVects[:,n_eigValIndice] #最大的n个特征值对应的特征向量 lowDDataMat=newData*n_eigVect #低维特征空间的数据 reconMat=(lowDDataMat*n_eigVect.T)+meanVal #重构数据 return lowDDataMat,reconMat #newData,meanVal=zeroMean(X) lowDDataMat,reconMat=pca(X,n=1) p=[] for i in range(24): j=i#np.log10(i) p.append(j) plt.plot(p, lowDDataMat, color='r', linewidth=0.2, alpha=0.6,label="after PCA") #plt.plot(p, X, color='b', linewidth=0.2, alpha=0.6,label="init") plt.legend()#显示图例名字 #pdf保存 plt.savefig('1.pdf') plt.show()

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162