有个PCA的东西-白红宇

有个PCA的东西

阅读量：4148 次

发布时间：2019-05-25

本文共 8933 字，大约阅读时间需要 29 分钟。

一、PCA入门

最近在看pca，pca可以用于将维，下面通过举一个具体的例子还说明pca的使用，为了便于可视化，我们假设数据是2维的。数据来源：

上面的是一个2维的数据x=(x1,x2)，一个45对。下面绘制出这个45对数据在的图如下：

图1

第一步：对数据进行预处理，使得特征x1和x2的相应的均值为0，即数据中心化，然后通过公式：

，

计算协方差。

第二步：求解协方差矩阵的特征值和特征向量，公式为A*x=lambda*x，这里面可以通过svd求解得到特征向量，通过计算得到特征向量为U=[u1,u2]，其中

u1=[-0.7071,-0.7071]'，u2=[-0.7071,0.7071]'，对应的特征值为：0.1607和0.0153。这里面u1表示主要特征向量，u2表示第二个特征向量。由于协方差矩阵是一个实对称矩阵，所以它的不同特征值对应的特征向量是正交的，因此在这里u1和u2是正交的，即u1和u2的点积为0，同时u1和u2的模都为1，即|u1|=|u2|=1。下面将向量u1和u2也绘制到图上。图如下：

图2

第三步：对于每个样本x，，表示x在向量u1上的投影长度，这是因为u1x=|u1||x|cos(theta)，故x在u1上的投影为：u1x/|u1|=|x|cos(theta)，因为|u1|=1，因此u1x=u1'x=|x|cos(theta)。同理x在u2上的投影长度为：u2'x，由于u1和u2是是正交的，可以以(u1,u2）为基底表示x，这样有x=(u1'x)u1 + (u2'x)u2成立，很显然，在以u1,u2为基底的坐标系中x可以使用[u1'x,u2'x]'来表示，下面是对于45对数据的在以u1,u2为基底的坐标系上的投影，图如下：

图3

在上图中，坐标系上的每个点中，对于某一个x而言，图中的横坐标是是u1'x的值，图中的纵坐标是u2'x的值，我们可以使用如下公式表示：

$\begin{align}x_{\rm rot} = U^Tx = \begin{bmatrix} u_1^Tx \\ u_2^Tx \end{bmatrix} \end{align}$ ，

这样对于所有的训练样本有： $\textstyle x_{\rm rot}^{(i)} = U^Tx^{(i)}$ ，所有可以得到上面的图。

由于 $\textstyle U$ 是一个正交矩阵，即 $\textstyle U^TU = UU^T = I$ ，所以可以从 $\textstyle x_{\rm rot}$ 到 $\textstyle x$ ，进行转换，公式为： $\textstyle U x_{\rm rot} = UU^T x = x$ 。

第四步：数据降维：由图1和图3可以看出我们给定的45对数据的主要方向是u1或 $\textstyle x_{{\rm rot},1}$ 是的方向上，从图3我们发现 $\textstyle x_{\rm rot}$ 分布在x_rot,2=0轴附近，这样如果我们需要

把 $\textstyle x_{\rm rot}$ 降低到1维的，那么我们应该有：

$\begin{align}\tilde{x}^{(i)} = x_{{\rm rot},1}^{(i)} = u_1^Tx^{(i)} \in \Re.\end{align}$

更一般的，如果 $\textstyle x \in \Re^n$ ，如果我想把x降到k维，那么我们可以取 $\textstyle x_{\rm rot}$ 得k个主件。如果我们考虑 $\textstyle \tilde{x}$ 和 $\textstyle x_{{\rm rot}}$ 的维度一样，那么应该可以写成这样：

$\begin{align}\tilde{x} = \begin{bmatrix} x_{{\rm rot},1} \\\vdots \\ x_{{\rm rot},k} \\0 \\ \vdots \\ 0 \\ \end{bmatrix}\approx \begin{bmatrix} x_{{\rm rot},1} \\\vdots \\ x_{{\rm rot},k} \\x_{{\rm rot},k+1} \\\vdots \\ x_{{\rm rot},n} \end{bmatrix}= x_{\rm rot} \end{align}$ ，

其实上面的公式中，本质上就是使得 $\textstyle U$ 的第k+1项后的所有向量为0向量。下面是 $\textstyle \tilde{x}$ （ $\textstyle n=2, k=1$ )的图像

图4

从图3到图4，很显然，这是相当于把图3所有坐标点的纵坐标变为了0，从而得到了图4。

第五步：恢复估计的数据：有前面可知 $\textstyle x = U x_{\rm rot}$ ，那么由 $\textstyle \tilde{x}$ 如何恢复到 $\textstyle \hat{x}$ ，显然我们可以使用如下公式：

$\begin{align}\hat{x} = U \begin{bmatrix} \tilde{x}_1 \\ \vdots \\ \tilde{x}_k \\ 0 \\ \vdots \\ 0 \end{bmatrix} = \sum_{i=1}^k u_i \tilde{x}_i.\end{align}$

下面是 $\textstyle \hat{x}$ 的图：

图5

我们可以对比一下图1和图5，很明显数据都被投影到了u1上了。

第六步：对于多维的情况，降到多少维(k)比较合理呢？由给定数据的协方差 $\textstyle \Sigma$ 可以得到它的特征向量 $\begin{align}U = \begin{bmatrix} | & | & & | \\u_1 & u_2 & \cdots & u_n \\| & | & & | \end{bmatrix} \end{align}$ ，以及它对应的特征值 $\textstyle \lambda_1, \lambda_2, \ldots, \lambda_n$ 。当 $\textstyle k=n$ ，可以100%的估计数据，当 $\textstyle k=0$ ，可以0%的估计，为此为了确定k，我们可以定义： $\begin{align}\frac{\sum_{j=1}^k \lambda_j}{\sum_{j=1}^n \lambda_j}.\end{align}$ ，对于2D的数据，如上面lambda1=0.1607和lambda2=0.015，那么对于k=1，0.1607/(0.1607+0.015)= 0.9146，通常情况下，当 $\begin{align}\frac{\sum_{j=1}^k \lambda_j}{\sum_{j=1}^n \lambda_j}.\end{align}$ 的值大于85%也OK了。

代码实现：

close all    %%================================================================  %% Step 0: Load data  %  We have provided the code to load data from pcaData.txt into x.  %  x is a 2 * 45 matrix, where the kth column x(:,k) corresponds to  %  the kth data point.Here we provide the code to load natural image data into x.  %  You do not need to change the code below.  %x的一列代表一个样本，一个样本有2个属性，即是2维数据  x = load('pcaData.txt','-ascii');  figure(1);  scatter(x(1, :), x(2, :));  title('Raw data');  xlabel('x1');  ylabel('x2')    %%================================================================  %% Step 1a: Implement PCA to obtain U   %  Implement PCA to obtain the rotation matrix U, which is the eigenbasis  %  sigma.     % -------------------- YOUR CODE HERE --------------------   u = zeros(size(x, 1)); % You need to compute this  %u=[u1,u2]  [n,m]=size(x); %n:维数;m：样本个数  x = x - repmat(mean(x,2),1,m); %样本中心化,均值为0  %其实m除和不除都没关系，最终计算出来的特征向量都是一样的，唯一的区别是特征值不一样%特征值之间差别也是差了1/m而已sigma = (1.0/m)*x*x';          %协方差  cov   = (1.0/m)*x'*x;%sigma和cov有一部分特征值是相同的[aa,bb,cc]=svd(cov);%u是特征向量，s是特征值[u,s,v] = svd(sigma);  %a是特征向量，b是特征值  [a,b]=eig(sigma);      % ------------绘制特征向量[u1,u2]------------------------------------   hold on  plot([0 u(1,1)], [0 u(2,1)]);  plot([0 u(1,2)], [0 u(2,2)]);  grid  scatter(x(1, :), x(2, :));  hold off    %%================================================================  %% Step 1b: Compute xRot, the projection on to the eigenbasis  %  Now, compute xRot by projecting the data on to the basis defined  %  by U. Visualize the points by performing a scatter plot.    % -------------------- YOUR CODE HERE --------------------  %计算所以点在以特征值为基上的投影，也就是x在[u1,u2]空间中的位置xRot = zeros(size(x)); % You need to compute this  %需要注意的是,|u1|和|u2|都为1，u=[u1,u2]%以u作为特征空间xRot = u'*x; %原始数据在u上的投影  %另外x可以从xRot中恢复：xRot = u'*x -> u*xRot = u*u'*x = (u*u')*x=x  % -------------------绘制旋转后的x值（此时x处于u1,u2的空间）-----------------     % Visualise the covariance matrix. You should see a line across the  % diagonal against a blue background.  figure(2);  scatter(xRot(1, :), xRot(2, :));  xlabel('x_r_o_t_,_1')  ylabel('x_r_o_t_,_2')  grid  title('xRot');    %%================================================================  %% Step 2: Reduce the number of dimensions from 2 to 1.   %  Compute xRot again (this time projecting to 1 dimension).  %  Then, compute xHat by projecting the xRot back onto the original axes   %  to see the effect of dimension reduction    % -------------------- YOUR CODE HERE --------------------   % 将数据进行降维，这里从2维降到1维k = 1; % Use k = 1 and project the data onto the first eigenbasis  xHat = zeros(size(x)); % You need to compute this  %xHat = u*([u(:,1),zeros(n,1)]'*x);  %这里是把纵坐标消除了。xrot =u'x =[u1'x;u2'x] --> [u1'x; 0],这是因为u2对应的特征值很小v2xrot1=([u(:,1),zeros(n,1)]'*x);  xrot2=[u(:,1)]'*x;%数据恢复,这里并没有添加0项，效果和xHat是一样的。%即xHat1 == xHatxHat1 =[u(:,1)]*xrot2;xHat = u*xrot1;  figure(3);  scatter(xrot1(1, :), xrot1(2, :));  title('xrot1');  grid    % --------------------------------------------------------   figure(4);  scatter(xHat(1, :), xHat(2, :));  title('xHat');  grid  hold on  plot([0 u(1,1)], [0 u(2,1)]);  plot([0 u(1,2)], [0 u(2,2)]);

二、PCA进阶

在对数据进行主成分分析（PCA）的时候，有的时候我们会遇到样本的维数>>样本的个数，如在人脸识别中，通常样本的维数达到了上万个，而样本个数可能是几百个。这样，如果我们需要计算这些数据的协方差，那么意味着协方差的维数将特别高(样本维数*样本维数)，从而导致计算的时间复杂度很高（在使用matlab计算的时候会出现内存不够），但是我们的确需要这些数据协方差相应的特征值和特征向量，那我们应该如何处理这个问题呢？有一种叫做SVD（奇异值分解）可以帮我们计算得到。

A. 奇异值

首先，我们先介绍奇异值分解，设A是M*N的实矩阵，那么存在M阶正交阵U，N阶正交阵V，使得

A = U*S*V‘

其中，S=diag(σi,σ2,……,σr)，σi>0 (i=1,…,r)，r=rank(A)，σi是A的奇异值，与此同时实矩阵A还有以下性质：

其中，U是A*A’的单位(因为U的每一列向量的模为1)特征向量矩阵，V为A‘*A的单位特征向量矩阵。

对于任意矩阵A，它的奇异值是AA’或者A‘A的非零特征值的开方。

下面是一些测试程序

clear,clcformat short g%压缩空格format compact%M < N 情况A = reshape([1:18],3,6);[U,S,V] = svd(A);index =find(S < 1e-10);S(index)=0;U,S,VA1 = U*S*V'[U1,S1] = eig(A*A');index =find(S1 < 1e-10);S1(index)=0;U1,S1[U2,S2] = eig(A'*A);index = find(S2 < 1e-10);S2(index)=0;U2,S2sqrtS1=sqrt(S1)sqrtS2=sqrt(S2)%%%M > N的情况clear,clcA = reshape([1:18],6,3);[U,S,V] = svd(A);index =find(S < 1e-10);S(index)=0;U,S,VA1 = U*S*V'[U1,S1] = eig(A*A');index =find(S1 < 1e-10);S1(index)=0;U1,S1[U2,S2] = eig(A'*A);index = find(S2 < 1e-10);S2(index)=0;U2,S2sqrtS1=sqrt(S1)sqrtS2=sqrt(S2)

B.奇异值和特征值的关系

当矩阵是共轭对称矩阵的时候，矩阵的特征值=矩阵的奇异值，其他情况不一定相同，如下：

T = 30    70   110   150    70   174   278   382   110   278   446   614   150   382   614   846%%%通过对T进行奇异值分解得到，其中u和v是正交阵，s是奇异值u =  -0.13472     -0.82574      0.52489      0.15649     -0.34076     -0.42882     -0.81649      0.18258     -0.54679    -0.031892     0.058312     -0.83463     -0.75283      0.36503      0.23329      0.49556s =    1491.7            0            0            0            0       4.2904            0            0            0            0            0            0            0            0            0            0v =  -0.13472     -0.82574      0.33649     -0.43217     -0.34076     -0.42882     -0.12654      0.82704     -0.54679    -0.031892     -0.75641     -0.35756     -0.75283      0.36503      0.54645    -0.037308%%%通过对T进行特征分解得到，其中u1是特征向量矩阵，s1是特征值矩阵u1 = -0.49633      0.23163     -0.82574      0.13472      0.48913     -0.67878     -0.42882      0.34076      0.51074      0.66268    -0.031892      0.54679     -0.50354     -0.21553      0.36503      0.75283s1 =        0            0            0            0            0            0            0            0            0            0       4.2904            0            0            0            0       1491.7%%%通过以上可以得出，但T是对称矩阵的时候，通过svd得到的奇异值和eig得到的特征值的确相等，%%%而正交阵u和特征向量矩阵u1只是符号和顺序改变了而已。

C.奇异值分解和PCA关系

首先，我们要明白我们的问题是什么，假设A是我们的中心化后（均值=0）数据，它是一个样本维数>>>样本个数的数据，即A是一个MxN的矩阵，M是样本个数，N是样本维数，其中，N>>M。现在我们希望求解X的协方差矩阵A*A'(这里A是NxM的矩阵，实际上是协方差是1/(m-1)*A*A'，但这并不影响结果)，但是A*A'不好求解它的特征值和特征向量，但是A'*A的特征值和特征向量好求，同时他们有一个关系，那就是他们的非零特征值相同，此外它们还有下面的关系：

即可以通过计算A'*A的特征向量vi和奇异值σi（这里和特征值的关系上面有），可以得到A*A'的特征向量ui，这样我们的问题就得到了解决了，即我们得到A*A'的协方差矩阵的特征向量了！

下面我们直接给出一个基于PCA的人脸识别程序：

%PCA_SVM人脸识别clear,clcglobal im;%使用全局变量imgdata=[];%训练图像矩阵N=15;  %人数M=1;  %每人多少张for i=1:N  for j=1:M    a=imread(strcat('.\orl_faces\s',num2str(i),'\',num2str(j),'.bmp'));    m=size(a,1);    n=size(a,2);    b=a(1:size(a,1)*size(a,2)); % b是行矢量 1×（图像row*图像col）    b=double(b);    imgdata=[imgdata; b]; % imgdata 是一个N * (图像row*图像col)   end;end;imgmean = mean(imgdata); %平均图片,1x(图像row*图像col)minus=imgdata-repmat(imgmean,N*M,1);% minus是一个N*(图像row*图像col)矩阵，是训练图和平均图之间的差值(样本中心化）%这个地方有一个小技巧，M'*M和M*M'的特征向量相同%sigma = minus'*minus;%[aa,bb,cc]= pcacov(sigma);covx=minus* minus'; % N * N (样本个数*样本个数)阶协方差矩阵 %PCA，用协方差矩阵的转置来计算以减小计算量, %coeff投影矩阵（由特征向量构成U=[u1,u2,..,un],并且已经按特征值大小排好了对）， %u1和u2,...un正交，且他们的模都为1 %latent特征值，explained：特征值的累积比例  [coeff, latent,explained] = pcacov(covx);% 选择构成95%的能量的特征值%选择特征向量i=1;proportion=0;while(proportion < 95)  proportion=proportion+explained(i);  i=i+1;end;p=i-1;% 训练得到特征脸坐标系i=1;%降维的过程%这里得到的base其实就是数据样本的协方差的特征向量矩阵(已经选出来了，即降维)while (i<=p && latent(i)>0) % base是N×p阶矩阵，用来进行投影，除以latent(i)^(1/2)是对人脸图像的标准化， % covx1=minus*minus'简化后的特征向量  base(:,i) = latent(i)^(-1/2)*minus' * coeff(:,i);   i = i + 1;endreference = imgdata*base; %将图像数据投影到base的空间中，得到降维的数据xRot = base'*imgdata';X=imgdata';%特征量for i=1:12    img2 = reshape(base(:,i),m,n);    figure(i);    filename = sprintf('特征脸%d.jpg',i);    imwrite(img2,filename);    imshow(img2,[]);end%数据恢复,即X是恢复后的数据（但是已经压缩了）XX=base*xRot;for i=1:15img1=reshape(XX(:,i),m,n);%figure(i)filename = sprintf('压缩后的脸%d.jpg',i);imwrite(img1,filename);%imshow(img1)end

原始图像数据